r 


AD-757  880 


ORBITAL  PARAMETER  DETERMINATION  BY 
WEIGHTED  LEAST  SQUARE  ERROR  AND 
KALMAN  FILTERING  METHODS 

Joseph  J.  Pollard 

Air  Force  Institute  of  Technology 
W r i g h t - P a t t e rs o n  Air  Force  Base,  Ohio 

December  1972 


J 


r 


DISTRIBUTED  BY: 


National  Technical  Information  Service 
U.  S.  DEPARTMENT  OF  COMMERCE 

5285  Port  Royal  Road,  Springfield  Va.  22151 


THIS  DOCUMENT  IS  BEST 
QUALITY  AVAILABLE.  THE  COPY 
FURNISHED  TO  DTIC  CONTAINED 
A  SIGNIFICANT  NUMBER  OF 
PAGES  WHICH  DO  NOT 
REPRODUCE  LEGIBLY. 


Reproduced  by 

NATIONAL  TECHNICAL 
INFORMATION  SERVICE 

U  S  Department  of  Commerce 
_ Springfield  VA  22151 


•  t  ■ 

-  .<*  .  » 

*  ’  "  #.‘9  ► 

•  '  V'  •  '*•  V 

-  vV-. 

.  ■*  .. 

•V- 

.  •  >1 

.  VtT- 

f.  .!•  ....  :  « 
'  '  1 

•  y; ' , 

■  *-  i 

>.i,  >:•  .■ 

•  ...  s 

.  \} 

■  ■:  •>' 

■  v.*» 

*•’  A;. . 

.  ‘  .  •  -  /'*.  •  • 

r  -■  v 

- S 

w  ■ 

-  V  VJt.  ,  .  . 

■k.Hr«  • 

'*  •  '  .  ’•*  f 

UNCLASSIFIED 

Security  Classification 


DOCUMENT  CONTROL  DATA  -  R  &  D 

(Security  ciassillcation  ol  titie,  body  ol  abstract  and  indexing  annotation  must  be  entered  when  the  overall  report  Is  classllled) 


t.  ORIGINATING  ACTIVITY  (Corporate  author ) 


Air  Force  Institute  of  Technology  (AFIT/ENE) 
Wright-Patterson  Air  Force  Base,  Ohio  45433 


2*5.  REPORT  SECURITY  CLASSIFI^AT  DN 

Unclassified 


2b.  GROUP 


3.  REPORT  TITLE 


Orbital  Parameter  Determination  by  Weighted  Least  Square  Error  and  Kalman  Filtering 
f  Methods 

O'  i  j l  . j 1  .  ..;i  ' .  • 


4.  DESCRIPTIVE  NOTES  (7 "ype  ol  report  and  Inclusive  dates) 

AFIT  Thesis 

■::ir;  .  .  i  !  . 

.  J'  .1 

8  AUTHOR(S)  (First  name,  middle  initial,  last  name) 

Joseph  J.  Pollard 

Captain  USAF 

.  .  1  .  .  *  \A  .  i  .  . . 

.  .  ”  J 

8.  REPORT  DATE 

7a.  TOTAL  NO.  OF  PAGES 

7b.  NO.  OF  REFS 

December  1972 

145 

•  .  21  •  ‘ 

8a.  CONTRACT  OR  GRANT  NO. 

Da.  ORIGINATOR’S  REPORT  NUMBER(S) 

b.  PROJEC  T  NO. 

GGC/EE/73-13 

I  c  N/A 

j  d- 

Ob.  OTHER  REPORT  NO(S)  (Any  other  numbers  that  may  be  assigned 
this  report) 

I  10.  OISTRI  QUTION  STATEMENT 

This  document  has  been  approved  for  public  release  and  sale;  its  distribution  is 
limited. 

tt.  SUPPLEMENTARY  NOTES 

\2.  SPONSORING  MILITARY  ACTIVITY 

13.  ABSTRACT 


The  Extended  Kalman  Filter  and  Weighted  Least  Square  Error  filtering  techniques 
are  used  to  determine  estimates  of  the  classical  orbital  parameters  of  a  passive 
near-earth  satellite.  Modelling  of  the  geopotential  is  complete  through 
Modelling  of  the  earth's  atmosphere  ircludes  day-night  variations  as  well  as  altitude 
variations.  Kalman  smoothing  is  also  performed.  Both  filter  techniques  were  used 
on  a  simulated  orbit  as  well  as  a  polar  circular  orbit  and  a  low  perigee  eccentric 
orbit  obtained  from  actual  radar  data.  The  determination  of  the  orbital  parameters 
of  the  simulated  orbit  yielded  an  absolute  as  well  as  a  relative  comparison  of  the 
filter  techniques.  The  analyses  of  the  orbits  determined  from  actual  data  yielded 
a  relative  comparison  of  the  filter  techniques.  Comparable  results  were  achieved 
with  both  filter  techniques.  Analyses  include  discussion  of  orbital  perturbations 
as  well  as  mean  orbital  elements.  All  orbits  considered  were  multi-pass  with  a 
limited  number  of  observations  per  revolution. 


DD  ,Frj 473  \a  UNCLASSIFIED 

Security  Classification 


I  4. 


UNCLASSIFIED 

Security  Classification 


KEY  WORDS 


LINK  B 


-  ”  3 

WT 

ail 

ViT 

Orbit  Determination.., . .  . 

• . .  : 

0  t  f.  - 

V.- 

o  , 

.  J 

Atmospheric  Modelling  ■"-*  >• 

’/  ^  " ! 

Geopotential  Modelling  :n-T. 

i  :i 

Z  .  .  *! 

.*  j 

Extended  Kalman  Filter 

Weighted  Least  Square  Filtering 

Kalman  Smoothing 

i.i  I-./; 

Filtering 

i 

Orbital  Perturbations  /i.f 

_  t  # 

-r'.r  ...:\jnD 

/  \  . 

'  r  ;  \r;:T  t  ..  :  :  y<  .  in  *r 

i  \,*J  . 

• 

-• 

it  7,. . 

■;  o 1 1 1  f  ."cut:--!  \  ho  ♦  *  *  i  i  o  .*1  .... 

■  'j  O’* 

:  :  rti; 

■  ):Ij  S  •’ 

r- :>.?.  :■<  *  •  ij  urt»t:  '.  oor  • 

■>  < 

i  7 

‘V  .  . 

‘  i  J  •  j  • 

!.  r  .oo  ir.Iv::.;  ••  j  r  t 

ri . 

j  ; 

.  fj 

•’  •' 

ov  ir’In  -  >:  i  .  -j .  /  ;.t;  a.'  <•'  .  ;  1  s  .-.r.i.  c  >!. ..  i  on  i  „*• 

•J  ".‘i 

ji;  r. 

*  .  j 

. ; ;  '  0 

::oiipr  .r.  ■>-  ".oJu  r:r o  t  . u-j  t  ,'rrjn  o  :• 

.  .  1 1'j'j 

o  r-r :  :  » u  ;*'  f~u"  oi  ;  7t<  r  luo  -  .o 

I*, 

I  J.'.  .1 

.  Jr  ;; 

;  .J.'j. 

vfjt  .  "f  :•  fan  .  -if  ;:o?  : /.  •:■.  !'  ..••Jr.i 

*»  * .  #'  • 

*  J 

m  .*  or, :;•(.;•*  :c»j  ov‘ i:  I  j*f  o  yr.  [  I  r.r  -r.  ?i; I •  ^ >  ;■ 

iff  ..o' 

1  io.'l 

.  ;..li 

• !o  I  r  !■?■.'.  ir.nl’it:  *n‘‘  Do”. '  n.Joi.  -  *..*  t 

*:0 

;  .  •  V. 

i. :  1  . . 

•  .■  ■  ,  ^‘v.  ■  :■>.)  ./v.i 

.  i ;  ■ 

U’  .'. ' 

i  ...  i 

•ti  1 .  T".  '■[•"!  . '  ••  :  .  -  T  f  "*0  .  i  ;.:U:.  •  oi.ulo;:  K  lit'/ 

mi/. 

i 

v  .  . 

V  ,  f  . 

Ji  ••.?’.  i:  /f ;  ..joc.'.-'...  0  .  .. 

. 

. :  .  . . 

j  i  vi ' 

::  ; 

.  ;r.-  r  t  = :  I  .-737  • 

0" 

i  0  ,  ’1 

.  *-'io  : 

r'A.Ij.'l 


/b 


unclassified 

Security  Classification 


I 


ORBITAL  PARAMETER  DETERMINATION  BY 
WEIGHTED  LEAST  SQUARE  ERROR  AND 
KALMAN  FILTERING  METHODS 
THESIS 


GGC/EE/73-13 


Joseph  J.  Pollard 
Captain  USAF 


m.  a  j| 

v  e 


ciaafiiiM 

'S  1S75 

iBSErire, 


i 


This  document  has  been  approved  for  public  release 
and  sale;  its  distribution  is  unlimited. 


& 


K 


ORBITAL  PARAMETER  DETERMINATION  BY 
WEIGHTED  LEAST  SQUARE  ERROR  AND 
KALMAN  FILTERING  METHODS 

THESIS 

Presented  to  the  Faculty  of  the  School  of  Engineering 
of  the  Air  Force  Institute  of  Technology 

Air  University 

in  Partial  Fulfillment  of  the 
Requirements  for  the  Degree  of 

Master  of  Science 
by 

Joseph  J.  Pollard,  B.S.E.E. 

Captain  USAF 

Graduate  Guidance  and  Control 

December  1972 


This  document  has  been  approved  for  public  release 
and  sale;  its  distribution  is  unlimited. 


GGC/EE/73-13 


Preface 

This  thesis  is  a  continuation  of  the  work  done  by  other  students  at 
the  Air  Force  Institute  of  Technology.  The  most  recent  of  these  theses 
was  by  Captain  Jackson  R.  Ferguson,  Jr.,  GA-72.  A  sincere  respect  for 
their  previous  work  has  been  cultivated  over  the  past  year,  and  a  sincere 
debt  of  gratitude  is  owed  them. 

I  would  also  like  to  express  my  appreciation  and  thanks  to  Lt.  Col. 
Russell  A.  Hannen,  AFIT-ENE,  my  advisor,  for  his  advice  and  instruction 
in  stochastic  control  theory.  I  also  gratefully  acknowledge  the  assis¬ 
tance  of  Captain  Jackson  R.  Ferguson,  Jr.,  1st  Aerospace  Con  Squadron, 

ENT  AFB,  Colorado  in  obtaining  actual  radar  data  and  analysis  of  it. 

Last  but  far  from  least,  I  wish  to  thank  my  wife.  Hike,  for  her  help  and 
deep  understanding  over  the  year  during  which  this  thesis  was  prepared. 

Joseph  J.  Pollard 


GGC/EE/73-13 


Contents 


Page 

Preface  .........  .  .  .  ii 

List  of  Figures . . . v 

List  of  Tables . viii 

List  of  Symbols . .  . . i* 

Abstract . . . 

I.  Introduction.  .........  .  1 

Statement  of  the  Problem . „ . 2 

Assumptions  and  Limitations  .  2 

Plan  of  Development . . . .  2 

II.  Modelling  and  Dynamics . . . 4 

Earth's  Physical  Shape . 4 

Earth's  Geomagnetic  Field  .  5 

Earth's  Atmospheric  Density  .  8 

Geometric  Differentiation  .  6 

The  Satellite . 9 

The  Radar . . . 

Coordinate  Systems . 

Dynamics  -  State  Equations . “ 

Initial  Conditions . *2 

Measurements . *4 

III.  Kalman  Filtering  16 

Prediction . . . 

Filtering . *8 

Smoothing . *9 

Reverse  Integration  .  29 

Class  I  Filter  Modification . .  21 

IV.  Weighted  Least  Square  Error  Filtering  .  22 

V.  Simulation  Study . 2^ 

The  Simulation  Orbit . .  2^ 

Simulation  Data  Generation.  .  . . 29 

Orbital  Perturbations  . . 39 

Secular  Perturbations . 39 

Oblate  Earth  Perturbations . 32 

Filter  Operation . 38 

Kalman  Smoothing . D* 


•  •  • 
in 


GGC/EE/73-13 


Page 


Fixed  Interval  Smoothing  61 

Fixed  Point  Smoothing . 61 

Covariance  Matrix  Behavior  .  63 

Initial  Covariance  Matrix,  P(0) . 63 

Experimentation  with  the  Kalman  Gain,  K .  68 

Experimentation  with  the  value  of  P . 68 

Experimentation  with  the  value  of  Q . 69 

Pointing  Information  from  Prediction  .  69 

VI.  Actual  Data:  Ent  Set  I . 71 


Radar  Characteristics . 

Comparison  of  Orbital  Parameters 

Kalman  Smoothing  . 

Interpass  Pointing  Data . 

Filter  Covariance . 

Filter  Effectiveness  . 

VII.  Actual  Data:  Ent  Set  II  . 

Data . 

Radar  Characteristics . 

Comparison  of  Orbital  Parameters 

Kalman  Smoothing  .  . 

Covariance  of  the  Filters.  .  .  . 
Interpass  Pointing  Data . 

VIII.  Conclusions  and  Recommendations.  . 

Simulation  Conclusions  .  .  .  .  . 

Actual  Data  Conclusions . 

Recommendations . 


Bibliography  .  123 

Appendix  A:  Derivation  of  the  Gravitational  Field  of  the  Earth.  .  125 

Appendix  B:  Direction  Cosine  Matrices  .  128 

Appendix  C:  Derivation  of  the  System  Matrix,  F .  129 

Appendix  D:  Measurement  Sensitivity  Matrix,  M  .  131 

Appendix  E:  Gaussian  Distribution  Density  Curves . 132 

Appendix  F:  Determination  of  0_ . 141 

go 

Appendix  G:  Determination  of  Classical  Orbital  Elements  from  the  142 
State  Vector . . . 

Vita:  Joseph  J.  Pollard . 145 


73 

73 

81 

84 

84 

84 

90 

90 

90 

90 

110 

110 

117 

119 

119 

120 
121 


GGC/EE/73-13 


List  of  Figures 


Figure  Page 

1  Density  Modelling . , . 8 

2  Earth  Centered  Inertial  and  Rotating  and  Station  Centered  13 

Topocentric  Coordinate  Systems  .  .  „  . . 

3  Gaussian  Density  Distribution  500  Points  .  27 

4  True  Simulated  Semi-Major  Axis . 31 

5  True  Simulated  Eccentricity . 32 

6  True  Simulated  Angle  of  Inclination . ,  ,  .  .  .  33 

7  True  Simulated  Line  of  Nodes . 34 

8  True  Simulated  Period . 35 

9  Semi-Major  Axis  Error  Pass  1 . 42 

10  Semi-Major' Axis  Errors  Pass  2 . 43 

11  Eccentricity  Errors  Pass  1 . 44 

12  Eccentricity  Errors  Pass  2 . 45 

13  Angle  of  Inclination  Errors  Pass  1 . 46 

14  Angle  of  Inclination  Errors  Pass  2 . 47 

15  Line  of  Nodes  Errors  Pass  1 . 48 

16  Line  of  Nodes  Errors  Pass  2 . 49 

17  Angle  at  Perigee  Errors  Pass  1  .  . . 50 

18  Angle  at  Perigee  Errors  Pass  2 . 51 

19  Height  at  Perigee  Errors  Pass  1 . *  52 

20  Height  at  Perigee  Errors  Pass  2 . 53 

21  Period  Errors  Pass  1 . 54 

22  Period  Errors  Pass  2 . .  .  55 

23  Kalman  Determined  Semi-Major  Axis,  Simulated  Orbit  ....  56 

24  Kalman  Determined  Eccentricity,  Simulated  Orbit . 57 


v 


GGC/EE/73-13 


Figure  Page 

25  Kalman  Determined  Angle  of  Inclination,  Simulated  Orbit.  .  58 

26  Kalman  Determined  Line  of  Nodes,  Simulated  Orbit  .  59 

27  Kalman  Determined  Period,  Simulated  Orbit . 60 

28  Covariance  Simulated  Data  ?xx.  Simulated  Orbit  .  64 

29  Covariance  Simulated  Data  P^,  Simulated  Orbit  ......  65 

30  Standard  Deviation  of  Position,  Simulated  Orbit . 66 

31  Standard  Deviation  of  Velocity,  Simulated  Orbit.  .....  67 

32  Semi -Major  Axis  Ent  Set  1 . .  76 

33  Eccentricity  Ent  Set  1 . 77 

34  Angle  of  Inclination  Ent  Set  1 . .  78 

35  Line  of  Nodes  Ent  Set  1.  .  .  . . .  79 

36  Period  Ent  Set  1 . 80 

37  Covariance  of  P  ,  Ent  Set  1 . 86 

38  Covariance  of  P^,  Ent  Set  1 . 87 

39  Standard  Deviation  of  Position,  Ent  Set  1 . 88 

40  Standard  Deviation  of  Velocity,  Ent  Set  1 . .  .  89 

41  Semi-Major  Axis,  Kalman,  Ent  Set  2 . 94 

42  Eccentricity,  Kalman,  Ent  Set  2.  .  . . 95 

43  Angle  of  Inclination,  Kalman,  Ent  Set  2 . .  96 

44  Line  of  Nodes,  Kalman,  Ent  Set  2 . 97 

45  Period,  Kalman,  Ent  Set  2 . 98 

46  Semi -Major  Axis,  WLS,  Ent  Set  2 . 99 

47  Eccentricity,  WLS,  Ent  Set  2 . .  100 

48  Angle  of  Inclination,  WLS,  Ent  Set  2 .  101 

49  Line  of  Nodes,  WLS,  Ent  Set  2 . .  .  102 

50  Period,  WLS,  Ent  Set  2 .  1°3 


GGC/EE/73-13 


Figure  Page 

51  Semi-Major  Axis,  Kalman  Revised,  Ent  Set  2, .  105 

52  Eccentricity,  Kalman  Revised,  Ent  Set  2 .  106 

53  Angle  of  Inclination,  Kalman  Revised,  Ent  Set  2 .  107 

54  Line  of  Nodes,  Kalman  Revised,  Ent  Set  2 .  108 

55  Period,  Kalman  Revised,  Ent  Set  2 .  109 

56  Covariance  of  P  ,  Kalman  Revised,  Ent  Set  2 .  113 


57  Covariance  of  Pj^,  Kalman  Revised,  Ent  Set  2- .  114 

58  Standard  Deviation  of  Position,  Kalman  Revised,  Ent  Set  2  115 

59  Standard  Deviation  of  Velocity,  Kalman  Revised,  Ent  Set  2  no 


60  Gaussian  Density  Distribution  100  Points .  133 

61  Gaussian  Density  Distribution  200  Points .  134 

62  Gaussian  Density  Distribution  300  Points .  135 

63  Gaussian  Density  Distribution  400  Points . 136 

64  Gaussian  Density  Distribution  600  Points .  137 

65  Gaussian  Density  Distribution  700  Points .  138 

66  Gaussian  Density  Distribution  800  Points .  139 

67  Gaussian  Density  Distribution  900  Points .  140 


vii 


9 


| 

0 

1 


•] 


1 

& 

4 

i 


1 
J§ 


GGC/EE/73-13 


List  of  Tables 

Table  Page 

I  Theoretical  Radar  Parameters .  26 

II  Gaussian  Distribution  Data  Analysis .  28 

III  Orbital  Perturbation  Analysis  .  .  .  . .  36 

IV  Kalman  Filter  Analysis . .  40 

V  WLS  Filter  Analysis  .  41 

VI  Kalman  Smoother  Analysis . .  62 

VII  .  Theoretical  Interpass  Pointing  Data  Kalman  Filter  .  .  .  .  70 

VIII  Orbit  Data  Ent  Set  1 .  72 

IX  Radar  Parameters  Ent  Set  I.  .  *  . .  74 

X  Ent  Data  Set  I  Results .  75 

XI  Smoothing  Data  at  First  Observations  Ent  Set  1 .  82 

XII  Interpass  Pointing  Data  Ent  Set  I .  85 

XIII  Orbit  Data  Ent  Set  2 .  91 

XIV  Sensor  #344  Ent  Set  2  Parameters .  92 

XV  Ent  Data  Set  2  General  Perturbations .  93 

XVI  Ent  Data  Set  2  Special  Perturbations .  104 

XVII  Smoothing  Data  at  First  Observations  Ent  Set  2,  Special  111 

Perturbations . . . 

XVIII  Final  Covariances  Ent  Data  Set  2 . . .  112 

XIX  Interpass  Pointing  Data  Ent  Set  2  .  .  . .  118 


viii 


GGC/EE/73-13 


List  of  Symbols 


A 

B 

Cd 

F 

G 

I 

J 

Jk,m 

M 

N 

P 

pin 

k 

R 

S 

U 

V 

X 

Xs’VZs 

a 


ad  *ad  *ad 
x  y  z 

ag  »*ig  »ag 

sx  y 
e 

h 

i 

i.j.k 


State  Sensitivity  Matrix,  Weighted  Least  Squares 
Ballistic  Coefficient,  mZ/kgm,  .015  for  this  study 
Drag  Coefficient,  2.2,  for  this  study 
State  Sensitivity  Matrix 
Gravitational  attraction  of  the  earth 
Identity  matrix 

Cost  Functional,  Weighted  Least  Squares 
Coefficient  used  in  earth  modelling 
Measurement  Sensitivity  Matrix 
Index  used  with  Kalman  smoothing 
Covariance  of  the  states;  period 
Associated  Legendre  Polynomial 
Covariance  of  the  measurement  noise 
Frontal  area  of  the  satellite 
Earth's  geopotential 
Velocity 
State  Vector 

Coordinates  of  the  sensor  of  interest 
azimuth  angle 

Components  of  drag  acceleration 
Components  of  gravitational  attraction 

eccentricity;  elevation  angle;  geocentric  rotating  frame 
Non-linear  measurement  matrix 
Angle  of  inclination;  inertial  frame 
indices 


GGC/EE/73-13 


l 


m 

q 

t 

w 


y 

y 

E 


P 


* 

Q 

w 


l|F|| 

at 

A-l 

x(k|j) 

bo 


length  of  the  vehicle 
mass  of  the  vehicle 

-12 

constant  used  with  Level  I  modification,  1  x  10  usually 

time;  tangent  frame 

earth  rotation  rate 

measurement  vector 

point  of  Ares,  Vernal  Equinox 

criterion  of  convergence.  Weighted  Least  Squares 

offset  angle  between  inertial  and  rotating  geocentric 

East  Longitude 

range,  atmospheric  density 

State  Transition  Matrix;  North  Latitude 

Line  of  Nodes 

Argument  of  Perigee 

Vector  Quantity 

Norm  of  F  where  F  is  a  vector 

Transpose  of  matrix  A 

Inverse  of  matrix  A 

Value  of  x  at  kth  iteration  given  information  through  the 
jth  iteration 

b  evaluated  at  the  time,  t0 


GGC/EE/73-13 


Abstract 

The  Extended  Kalman  Filter  and  Weighted  Least  Square  Error  filtering 
techniques  are  used  to  determine  estimates  of  the  classical  orbital  para¬ 
meters  of  a  passive  near-earth  satellite.  Modelling  of  the  geopotential 
is  complete  through  Modelling  of  the  earth's  atmosphere  includes 

day-night  variations  as  well  as  altitude  variations.  Kalman  smoothing 
is  also  performed.  Both  filter  techniques  were  used  on  a  simulated  orbit 
as  well  as  a  polar  circular  orbit  and  a  low  perigee  eccentric  orbit  ob¬ 
tained  from  actual  radar  data.  The  determination  of  the  orbital  para¬ 
meters  of  the  simulated  orbit  yielded  an  absolute  as  well  as  a  relative 
comparison  of  the  filter  techniques.  The  analyses  of  the  orbits  deter¬ 
mined  from  actual  data  yielded  a  relative  comparison  of  the  filter  tech¬ 
niques.  Comparable  rssults  were  achieved  with  both  filter  techniques. 
Analyses  include  discussion  of  orbital  perturbations  as  well  as  mean 
orbital  elements.  All  orbits  considered  were  multi-pass  with  a  limited 
number  of  observations  per  revolution. 
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ORBITAL  PARAMETER  DETERMINATION  BY 
WEIGHTED  LEAST  SQUARE  ERROR  AND 
KALMAN  FILTERING  METHODS 
I.  Introduction ' 

When  the  first  vehicles  were  orbited  in  the  late  1950's  and  early 
1960's,  a  demand  for  greater  accuracy  and  higher  speed  in  the  computation 
of  orbital  parameters  became  evident.  Both  engineers  and  scientists  com¬ 
bined  their  knowledge  of  the  earth's  non-uniform  gravitational  field, 
the  earth's  oblate  spheroidal  shape,  and  the  non-linear  time  varying  at¬ 
mosphere  of  the  earth  to  produce  highly  accurate  models.  Communications 
engineers  strained  at  the  problem  of  more  accuracy  from  their  radars  while 
the  age  of  computer  technology  paced  the  speed  which  could  be  achieved  in 
the  determination  problem. 

Previous  AFIT  theses  (Ref:  6,7,  $  17)  have  addressed  many  of  the 
problems  of  orbital  determination  in  both  simplified  simulations  and 
quasi-real  world  environments  using  both  weighted  least  square  error  and 
Kalman  filtering  methods.  Success  has  been  demonstrated  in  the  applica¬ 
tion  of  these  filtering  methods  on  both  simulation  and  actual  radar  data 
measurements . 

However,  as  yet,  all  of  these  concepts  have  not  been  gathered  to¬ 
gether  using  extremely  accurate  modeling  or  even  exercising  the  filtering 
and  smoothing  capabilities  of  stochastic  theory  to  the  extent  possible. 
The  purpose  of  this  study  is  to  apply  both  the  Weighted  Least  Square 
Error  and  Kalman  filtering  theory  to  the  problem  of  real-world  near- 
earth  orbital  determination  using  extremely  accurate  models  of  the  earth 
and  its  atmosphere.  Smoothing  algorithms  will  be  used  in  conjunction 
with  the  Kalman  filtering  algorithms  to  determine  the  classical  orbital 
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parameters.  Perturbations  occuring  in  the  elements  will  be  investigated 
and  compared  with  theoretical  predictions. 

Special  interest  is  shown  in  this  study  for  the  short  term,  fast 
determination  problem  and  the  capability  to  predict  vehicle  position  and 
velocity  on  future  passes  in  order  to  assist  radar  acquisition,  provide 
pointing  data,  and/or  aid  intercept  capability.  The  multi-pass  multi¬ 
radar  problem  is  studied  assuming  limited  information  availability  on  all 
passes. 

The  problem  will  be  studied  by  first  addressing  highly  accurate 
models  to  a  simulated  orbit,  two  pass,  single  radar  problem.  Having  ob¬ 
tained  a  basis  for  comparison  of  both  filter  and  smoothing  capabilities 
on  data  with  known  true  parameters,  the  study  will  shift  to  analyses  of 
actual  radar  data.  The  first  case  to  be  studied  will  be  a  three  pass,  two 
radar  problem  of  a  satellite  in  near  circular  polar  orbit.  The  second 
case  will  be  a  limited  data  problem  using  only  one  radar  providing  measure¬ 
ments  of  a  satellite  having  an  eccentric,  inclined,  low  perigee  orbit.  In 
both  cases,  parameters  obtained  will  be  compared  to  parameters  obtained 
through  use  of  current  operational  capabilities.  Use  will  be  made  of 
both  general  and  special  perturbation  theory. 

During  this  study  the  following  assumptions  are  made: 

1.  The  satellite  under  consideration  is  spherical,  cylindrical,  or 
of  other  common  shape  with  no  large  solar  panels  or  massive  extended 
antennae  such  that  its  ballistic  coefficient  is  accurately  modeled  as 
B».015  m^/kgm  (Ref.  12:12-17). 

2.  The  satellite  is  non-thrusting  and  therefore  incapable  of  man¬ 
euvering. 

3.  The  earth's  gravitational  field  is  accurately  represented  by 
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inclusion  of  zonal,  t;sserai,  and  sectoral  harmonics  through 

4.  The  earth's  atmosphere  is  non-rotating  with  respect  to  the  surface 
of  the  earth.  The  complete  atmospheric  model  is  described  in  Chapter  2. 

So  The  earth';,  atmosphere  has  negligible  effect  above  800  kilo¬ 
meters  . 

I 

6.  The  earth's  atmosphere  has  sufficient  density  below  120  kilo¬ 
meters  that  orbital  degradation  and  re-entry  have  begun. 

7.  All  orbits  considered  occurred  during  a  period  of  mean  solar 
activity. 

8.  Uncertainties  in  latitude,  longitude,  and  height  of  the  tracking 
sensors  are  assumed  negligible. 

9.  The  noise  in  the  tracking  sensor  data  is  assumed  to  be  a  zero- 
mean,  white,  Gaussian  sequence. 

This  study  is  divided  into  eight  chapters.  The  second  chapter  pre¬ 
sents  a  detailed  description  of  the  models  used,  state  equations,  and 
initial  conditions.  Chapter  3  presents  the  Kalman  prediction,  filtering, 
and  smoothing  algorithms  used.  Chapter  4  introduces  the  statistical  least 
square  error  filter.  The  fifth  chapter  presents  the  details  of  the  simu¬ 
lated  orbit  analysis  and  its  results.  Chapter  6  presents  an  orbital  anal¬ 
ysis  of  an  actual  vehicle  in  a  polar  circular  orbit.  Chapter  7  deals  with 
an  eccentric,  inclined,  low  perigee  vehicle.  The  final  chapter  presents 
recommendations  for  future  study  along  with  the  conclusions  of  this  study. 
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II.  Modelling  and  Dynamics 

In  oi’dcr  to  accurately  determine  the  orbital  parameters  of  a  satellite, 
the  environment  which  the  vehicle  experiences  must  be  accurately  described 
and  modeled. 


The  Earth's  Physical  Shape 

For  many  years  science  has  been  i  vare  that  the  actual  physical  shape 
of  the  earth  is  not  spherical.  Current  geophysical  theory  models  the 
earth  as  an  ellipsoid.  An  ellipsoid  is  a  surface  of  revolution  having  a 
smaller  polar  diameter  than  its  equitorial  diameter.  From  the  analysis 
of  orbital  data,  many  values  have  been  determined  for  the  ellipticity  of 
the  earth.  However,  in  1964,  a  value  of  1./298.30  was  adopted  as  an  inter¬ 
national  reference  standard  (Ref:  3:172).  This  value  is  used  throughout 
this  study.  The  international  reference  standard  for  the  equatorial  radius 
of  the  earth  is  also  used.  Thus,  the  following  parameters  are  employed 
for  the  physical  shape  of  the  earth: 

f  =  298.30  (1) 

re  =  6,378,165  meters  (Ref:  3:171)  (2) 


Hie  Earth's  Geomagnetic  Potential 


The  earth's  gravitational  potential,  U,  is  given  by 

^  Pi,m(sin<j>)  /  . 

I  -  (  CMco5(.XE) 


u .  tii 


►  gX  clV  JL  IcLLJLUJlc 

/  6  /  k 

1  ♦  l  [  l 
\  k=2'm=0 


(3) 


rJl 


+  Sk,msinCmXE^ 

where  P£  (sin  0)  is  an  associated  Legendre  polynomial,  r  is  the  radius 
to  the  point  of  interest,  is  the  east  longitude  and  0  is  the  north 

latitude  (Ref  3:  174-175).  The  coefficients  C.  and  S,  are  determined 

by  a  combination  of  orbital  analysis  and  geodetic  survey.  Another  form  of 
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equation  (3)  uses  the  coefficients  .  However,  since  J.  =  fg~  +  s2 

,  K,m  k,m 

and  from  trigonometry  it  can  be  shown  that 


A  cos  0  +  B  sin  0  *=  C  cos  (0  +  a) 


where 

C  -  i2 

and 

a  =  tan"^-  (B/A) 

there  is  actually  no  difference  between  these  forms.  The  m  are  often 
referred  to  as  zonal  constants  if  m  =  0,  sectoral  constants  if  k  =  m,  and 
tesseral  constants  otherwise.  The  constants  used  in  this  study  are  found 
in  Reference  8. 

The  associated  Legendre  polynomial  is  given  by 


P?!(x)  =  (l-x2)®/2  £L  P.(x)  (4) 

dxm  K 

where  P^  (x)  is  a  Legendre  polynomial.  The  Legendre  Polynomials  used  in 
this  study  are  given  in  Appendix  A  (Ref:  13:870). 

Since  from  basic  physics  it  is  shown  that 

G  «  VU  (5) 

it  is  necessary  to  take  the  first  partial  derivatives  of  potential  with 
respect  to  x,  y,  and  z  in  the  geocentric  rotating  coordinate  frame  in  order 
to  define  the  gravitational  force  field.  Additionally,  the  second  partials 
with  respect  to  all  combinations  of  x,  y,  and  z  will  be  required  in  order 
to  calculate  the  state  equations.  The  derivation  of  the  first  partials 
is  given  in  Appendix  A.  The  second  partials  are  determined  by  geometric 
approximation  which  will  be  discussed  later  in  this  chapter. 


The  Earth’s  Atmospheric  Density 

The  basis  for  the  atmospheric  density  model  used  in  this  study  are 
the  1965  Cospar  (Council  on  Space  and  Atmospheric  Research)  atmospheric 
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data  tables  (Ref:  1:175-186).  The  earth's  atmospheric  density  is  found 
to  be  a  function  of  three  basic  variables:  (1)  height  above  the  reference 
ellipsoid,  (2)  position  with  respect  to  the  sun  or  local  time,  and  (3) 
solar  activity.  The  first  two  of  these  basic  variables  are  fully  accounted 
for  by  the  model  used  in  this  study  while  solar  activity  is  a  completely 
random  variable  and  is  treated  by  the  assumption  made  in  chapter  1  that  a 
condition  of  mean  solar  activity  is  present  during  all  orbital  calcula¬ 
tions. 

The  raw  data  relating  altitude,  local  time,  and  density  was  first 
fit  with  a  curve  relating  density  and  altitude.  After  these  coefficients 
were  determined  using  a  least  square  error  fit,  the  resulting  coefficients 
of  a  seventh  order  polynomial  were  fitted  versus  local  time  data.  Altitude 
increments  were  20  kilometers  and  time  increments  were  two  hours.  The 
time  versus  fitted  altitude/density  coefficients  fit  was  done  using  a 
fourth  order  curve.  Thus  an  analytic  model  was  formed  from  the  raw  data 
tables.  Curves  of  all  possible  orders  were  fitted  and  the  curve  of  least 
order  which  generated  a  good  fit  was  chosen. 

To  verify  that  no  undesired  peaks  or  discontinuities  had  occurred  due 
to  the  curve  fitting,  listings  by  altitude  in  increments  of  1  kilometer 
and  by  time  in  increments  of  3  minutes  were  made.  No  irregularities  were 
detected.  Figure  1  shows  groupings  which  were  made  during  the  curve  fitting 
analysis.  All  intersections  were  overlapped  at  the  boundary  points  to 
assure  continuity  when  the  analytic  model  was  formed. 

Geometric  Differentiation 

The  calculation  of  the  various  state  equations  which  will  be  discussed 
later  requires  the  evaluation  of  first  partials  of  the  density  as  well  as 
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second  partials  of  the  gravitational  potential.  To  accomplish  this 
geometric  differentiation  was  performed: 


where  f(x)  is  the  evaluation  of  the  function  at  x 

and  Ax  °  1  meter 

The  choice  of  the  delta  was  done  after  experimentation  showed  that  no 
change  occurred  when  a  smaller  delta  was  chosen.  This  choice  of  delta  is 
equivalent  to  one  part  in  one  million  for  applications  in  this  study. 
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Figure  1.  Density  Modelling 
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The  Satellite 

The  assumptions  that  the  satellite  is  non-maneuvering  and  has  a 
ballistic  coefficient  of  .015  have  previously  been  stated.  A  slight 
justification  for  the  latter  assumption  is  now  presented. 

King-llele  states  (Ref  12:12-18)  that  for  a  satellite  in  an  orbit 
with  eccentricity  between  0.0  and  0.2  the  value  of  Cd  ranges  from  between 

2.1  to  2.2  for  a  sphere  to  2.1  to  2.25  for  a  cylinder  inclined  to  the  air¬ 
flow.  For  a  cylinder  tumbling  end  over  end  is  determined  to  be  about 
2.15.  For  cones  with  a  serai-apex  angle  of  15  to  20  degrees,  is  near 
2.10.  Therefore  it  appears  that  Cd  may  be  assumed  to  be  approximately 

2.2  with  only  small  error. 

The  average  mean  cross  sectional  area  of  a  cylindrical  satellite  may 
be  expressed  by: 

S  =  i,d(.813  +  .25  d/Jt)  (7) 

where 

S  =  the  mean  cross  sectional  area 
l  =  the  length  of  a  cylindrical  shape 
d  =  the  diameter  of  the  cylinder 
Since  the  drag  force  on  a  vehicle  is  given  by 

Fd  *  -  1/2  CdSpV2  (8) 

Where  p  =  the  atmospheric  density 

V  =■  the  vehicle's  velocity 
and  the  other  terms  are  as  defined  above. 

The  ballistic  coefficient,  B,  is  defined  to  be 

B  =  CdS/M  (9) 

where  M  is  the  mass  of  the  vehicle.  Then  the  acceleration  associated 
with  the  drag  is  given  by: 
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a-d  =  -  BpV2/2  (10) 

In  a  typical  case,  using  the  above  equations  and  the  dimensions  of 
the  Discoverer  series  of  satellites  (Ref  11:182)  the  ballistic  coefficient 
is  determined  to  be 

B  =  .02  m2/kgm  (11) 

2 

which  is  sufficiently  close  to  the  assumed  value  of  .015  m  /kgm. 

Since  the  models  of  the  gravitational  field  and  the  atmosphere  are  for 
the  earth  this  study  is  necessarily  restricted  to  satellites  which  are  in 
near-earth  orbits. 

The  Radar 

The  radars  used  in  this  study  are  assumed  to  have  a  priori  biases  and 
sigmas  associated  with  the  measurements  of  range,  range  rate,  azimuth  and 
elevation  which  they  provide.  Sampling  periods  of  6  and  10  seconds  are 
used  in  this  study  since  they  are  typical  of  operational  radars. 

The  radar  is  assumed  located  at  a  Class  A  survey  point  from  which  we 
know  its  position  exactly. 

Coordinate  Systems 

In  this  study  computation  and  measurements  are  done  in  three  common 
coordinate  systems. 

The  first  of  these  is  the  inertial  system,  designated  by  i,  which  has 
its  x-axis  directed  toward  the  Vernal  Equinox  (or  the  Point  of  Aries),  the 
z-axis  through  the  north  geographic  pole,  and  its  origin  located  at  the 
geocenter.  For  our  purposes,  this  is  a  non-rotating,  non-accelerating 
computation  frame  where  Newton's  laws  of  motion  are  directly  applicable. 
This  is  the  frame  in  which  all  orbital  elements  are  calculated. 

The  second  coordinate  system  to  be  considered  is  the  rotating 
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geocentric  system,  e.  This  system  has  its  origin  colocated  with  the  i 
frame  at  the  geocenter;  its  x-axis  passes  through  the  Greenwich  meridian 
in  the  equitorial  plane.  The  z-axis  of  this  system  passes  through  the 
geographic  north  pole.  The  system  is  called  rotating  since  the  x-axis 
always  passes  through  the  Greenwich  meridian  rotating  at  approximately 
15°  per  hour.  Integration  of  the  state  equations  is  performed  in  this 
system. 

The  last  coordinate  system  which  we  consider  is  the  tangent  coordinate 
system.  The  origin  of  this  system  is  coincident  with  some  point  on  the 
surface  of  the  earth.  In  our  study  this  point  will  always  be  the  physical 
location  of  some  radar  of  interest.  This  localized  coordinate  system  has 
its  z  axis  opposite  in  direction  to  the  gravity  vector;  the  x  axis  points 
toward  North.  The  y-axis  then  points  West  completing  the  right-handed 
orthogonal  triad.  Measurements  taken  by  the  radars  are  first  expressed 
in  this  coordinate  system. 

Figure  2  shows  these  various  coordinate  systems  and  their  relative 
orientations. 

Since  transformations  from  each  of  these  frames  to  the  others  is 
necessary.  Appendix  B  gives  the  direction  cosine  matrices  involved. 

Dynamics — State  Equations 

With  the  description  of  the  vehicle  and  its  environment  complete,  we 
now  progress  to  the  dynamics  and  state  equations  relevant  to  the  current 
problem. 

The  States 

The  states  which  are  chosen  to  be  representative  of  the  orbital 
determination  problem  are  the  three  components  of  the  position  vector. 
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x,  y,  and  z  and  the  three  components  of  the  velocity  vector  x,  y,  and  z. 
Written  in  the  inertial  frame  the  state  equations  are: 


where  a  ,  a  ,  a  are  the  components  of  the  total  derivative  of  eqn  (3) 
gx’  gy  gz 

and  a^,  are  the  components  of  eqn  (10)*  Expressed  in  the  ro¬ 

tating  geocentric  frame  these  state  equations  become: 


where  w  =  earth  rotation  rate  and  where  the  laws  of  Coriolis  have  been 


applied  in  transferring  from  the  inertial  frame  to  the  geocentric  frame. 


Initial  Conditions 

Initial  conditions  will  have  to  be  supplied  so  that  the  integration  of 
these  equations  can  be  begun  regardless  of  what  method  of  filtering  is 
used.  Since  the  vehicle  may  not  enter  into  its  predicted  nominal  trajec¬ 
tory,  the  first  measurement  set  must  be  used  to  generate  the  initial  con¬ 
ditions.  However,  the  conditions  for  application  of  the  Kalman  filter 
require  that  the  initial  state  be  uncorrelated  with  the  measurement  noise. 
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Figure  2.  Earth  Centered  Inertial  and  Rotating  and  Station 
Centered  Topocentric  Coordinate  Systems 


13 


GGC/EE/73-13 


To  accomplish  this  the  azimuth  and  elevation  angles  are  rounded  to  the 
nearest  .25  degree.  The  range  is  rounded  to  the  nearest  1000  meters 
while  the  range  rate  is  rounded  to  the  nearest  10  meters  per  second.  Al¬ 
though  this  method  may  produce  a  larger  error  than  necessary  initial ly^ 
the  conditions  for  filter  application  are  met.  For  application  to  the 
Weighted  Least  Squares  filter,  the  nominal  value  used  is  the  value  deter¬ 
mined  directly  from  the  first  set  of  measurement  data. 

To  complete  the  required  set  of  six  initial  conditions,  the  azimuth 
and  elevation  rates  are  formed  by  use  of  equation  (6)  with  the  incremental 
time  equal  to  the  sampling  rate  of  the  radar  in  use  and  f(x  +  Ax)  occurring 
at  the  second  measurement.  Thus  six  initial  conditions,  azimuth,  azimuth 
rate,  elevation,  elevation  rate,  range  and  range  rate  are  supplied  for 
initializing  the  integration  of  the  state  equations. 

The  Measurements 

The  measurements  are  related  to  the  state  variables  through  the  follow¬ 


ing  equations: 

p  -  ((xx  -  xs)2  +  (x2  -  ys)2  +  (X3  -  ZJ)1/2  (14) 

P  B  -  ((x.  -  x  )x  +  (x  -  y  )y  +  (x  -  z  )x  )  (15) 

p  i  S  4  2  S  4  3  So 

a  a  tan’1  (-yt/xt)  (16) 

e  =  tan"1  (x  /(x£  +  y^) 1/2^  (17) 

l  W 


where  x^,  x2,  Xj,  x^,  x^,  x^=  the  vehicle  position  and  velocity 

expressed  in  the  rotating  geocentric 
frame,  the  state  vector 

x  ,  y  ,  z  ’  e  the  coordinates  of  the  tracking 

O  w  w 

sensor  in  the  rotating  geocentric 
frame 
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xt,  yt,  z£  =  the  coordinates  of  t.he  vehicle 

expressed  in  the  tangent  frame 
by  use  of  the  direction  cosine 
matrix  found  in  Appendix  B. 
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III.  Kalman  Filtering 

A  Kalman  Filter  is  an  algorithm  for  the  determination  of  the  states 
or  parameters  of  a  system  using  measurements  which  may  not  reflect  all 
of  the  states  or  parameters  directly  and  are  contaminated  by  white, 
Gaussian  noise.  This  rather  concise  statement  sums  up  the  intent  and 
purpose  of  this  chapter.  Equations  will  be  presented  by  which  the  Kalman 
filtering  algorithm  can  be  applied  to  the  orbital  determination  problem. 
No  rigorous  proof  of  the  equations  will  be  presented  since  they  have  al¬ 
ready  been  derived  and  rederived  numerous  times  (Refs:  10:Ch  8  and  14: 

Ch  4 ,  5 ,  and  6)  . 

All  discussions  of  Kalman  filtering  must  begin  with  a  discussion  of 
the  general  state  equations  descriptive  of  the  system: 

*  ■  f  (x,t,u,w)  (18) 

z  =  h  (x,t,v)  (19) 

where  x  is  the  state  vector 

u  is  a  deterministic  control  vector 
w  is  a  random  noise  vector  and 
v  is  a  measurement  noise  vector. 


These  are  non-linear  equations  which  in  general  are  extremely  cumbersome 
and  difficult  to  handle.  As  a  result  a  linearization  about  some  nominal 
path  or  trajectory  is  often  made  leading  to  equations  of  the  form: 


x  *  F  x(t) 

(20) 

z.(R) 

=  M  x(k)  +  v(k) 

(21) 

where 

and 

F  *  1C&L 
9x 

x  =  x(k  +  l|  k) 

(22) 

M(x)  _  3hW 

"  3£  x  =  x(k+l|k) 

(23) 
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The  exact  form  of  F  and  M  pertinent  to  the  orbital  determination  prob¬ 
lem  are  derived  in  Appendixes  C  and  D.  The  matrix  F  is  often  referred  to 
as  the  state  sensitivity  matrix  and  the  matrix  M  is  called  the  measurement 
matrix.  Since  a  linearization  was  made  to  obtain  equations  (20)  and  (21) 
we  must  maintain  small  perturbations  throughout  the  analysis.  To  do  this, 
the  present  state  is  always  determined  by  the  forward  integration  of 
equation  (13).  Thus  as  the  recursive  filter  begins  to  take  effect  all  of 
the  linearization  assumptions  used  become  more  valid. 

Prediction 

Prediction  very  simply  consists  of  integrating  the  state  equation 
forward  from  the  present  state  to  the  next  state  to  be  considered  through 
use  of  the  state  equations. 

t+At 

x(k  +  1 1 k)  =  x(k | k)  +  f  f(x(k|k))  dt  (24) 

t 

o 

where  k  =-1,0, 1,2 . 

The  state  error  covariance  which  is  defined  by 

T 

P(k  +  1 | k)  =  E(  (x(k  +  l|k)-x(k  +  1 | k) )  (x(k  +  l|k)-x(k  +  l|k)))(25) 
can  be  projected  forward  through  use  of 

P(k  +  1 | k)  =  0(k  +  1 | k)  P(k|k)  0T(k  +  1 | k)  (26) 

Where  0(k  +  l,k)  is  the  State  Transition  Matrix: 

00c  ♦  l.k)  ■  ^  11  (27) 

0(k  +  l,k)  is  determined  by  integrating  the  linear  differential  matrix 
equation 

0  (k  ♦  l*k)  =  F(x(k|k))  0(k  +  l,k)  (28) 

with  the  initial  conditions 

0(k,k)  -  I  (29) 
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Filtering.  Filtering  is  the  process  of  updating  the  best  estimate 
of  the  present  state  through  use  of  an  optimal  gain  matrix,  K,  and  the 
measurement  residuals,  Az.  The  measurements  are  predicted  by  applying 
the  value  of  the  present  best  estimate  of  the  states  to  the  measurement 
equation: 

£  =  h(x(k  +  1 | k))  (30) 

remembering  that  there  is  zero-mean  noise  in  the  measurements.  The 
measurement  residuals  are  then  formed  by  subtracting  the  best  estimate 
of  the  measurements  from  the  true  measurements. 

Az  =  £(k  +  1)  -  z(k  +  l|k)  (31) 

These  residuals  are  then  linearly  added  to  the  present  best  estimate  of 
the  state  after  multiplication  by  the  optimal  gain: 

x(k  +  l]k  +  1)  =  x(k  +  *ik)  +  K(K  =  1)  Az  (32) 

where  the  optimal  gain,  K,  has  been  determined  from 

K(k  +  1)  +  P(k  +  l|k)MT(MP(k  +  l|k)MT  +  R)"1  (33) 

We  observe  that  if  the  noise  covariance,  R,  is  large  the  resulting  effect 
is  to  diminish  K.  However,  if  R  is  small,  then  K  will  be  large.  This 
relation  could  be  anticipated  for  optimal  application  of  the  information 
contained  in  the  residuals.  The  error  covariance  of  the  corrected  state 


estimate  x(k  +  l|k  +  1)  is  determined  by  the  symmetric  equation: 
P(k  +  1 ] k  +  1)  =  (I  -  K(k  +  l)M)P(k  +  l|k)(I  -  K(k  +1)M)T  + 

K(k  +  1)  R  K(k  +  1)T 


(34) 


where  M  is  evaluated  at  time  t,  . . 

k+1 

The  filter  applied  in  this  manner  is  called  the  Extended  Kalman 
Filter  since  the  equations  are  linearized  about  the  current  state  estimate, 
x(k  +  l]k).  The  prediction  process  does  not  involve  any  linearity  assump¬ 
tions  since  the  non-linear  state  equations  and  the  non-linear  measurement 
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matrix,  M  are  used. 

Smoothing 

Smoothing  is  the  process  of  conveying  information  found  in  later 
measurements  back  to  the  values  of  states  which  have  been  previously 
determined.  There  are  three  basic  types  of  smoothing,  (1)  fixed  interval 
smoothing,  (2)  fixed  point  or  single  point  smoothing,  and  (3)  fixed  lag 
smoothing.  These  are  detailed  in  Ref  14,  chapter  6. 

'Two  types  of  smoothing  were  applied  to  the  orbital  determination 
problem  in  this  study,  fixed  interval  and  fixed  point  smoothing. 

To  perform  smoothing  requires  the  storage  of  many  quantities,  in¬ 
cluding  the  state  transition  matrix,  the  error  covariance  matrix  and  the 
states.  This  is  the  chief  disadvantage  from  a  computational  viewpoint. 

The  fixed  interval  smoother  does  its  work  after  the  problem  has 
been  finished  or  at  least  after  the  number  of  measurements  over  which 
smoothing  is  desired  have  been  taken.  The  smoother  is  basically  a 
’'reverse"  smoother  in  the  sense  that  it  takes  the  last  measurement  and 
projects  it  backwards  to  the  previous  measurement.  This  process  is  re¬ 
peated  until  the  origin  of  the  problem  is  reached  where  theoretically  all 
information  learned  during  the  process  has  been  projected  back  upon  the 
original  starting  states. 

The  optimal  fixed  interval  smoothed  estimate  x(k|N)  is  governed  by 
the  system  of  equations: 

x(k|N)  ■  x(k | k)  +  A(k)(x(K  +  l|N)  -  x(k  +  l|k))  (35) 

where  k  =  N-l,  N-2,....0  where  x(N|N)  is  the  boundary  condition  for  k=N-l. 
The  optimal  smoothing  gain  matrix  is  given  by 

A(k)  =  P(k|k)0T(k  +  l,k)P_1(k  +  ljk)  (36) 
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The  optimal  fixed  interval  smoothing  error  covariance  matrix  is  given  by 
P(k|N)  «  P(k[k)  ♦  A(k)(P(k  +  l[N)  -  P(k  +  l|k)AT(k)  (37) 

where  k  =  N-l,  N-2.....0  subject  to  the  boundary  condition  P(N/N)  for 
k»N-l. 

The  fixed  point  smoother  is  a  forward  smoother  which  after  it  has 
been  given  the  states  and  time  at  which  smoothing  is  desired  continually 
updates  the  values  of  the  states  and  covariance  errors  at  that  point  as 
new  information  becomes  available. 

The  system  of  equations  governing  the  fixed  point  smoother  are: 
x(k|j)  =  x(kjj-l)  ♦  B(j)(x(j[j)-x(j|j-l)) 
for  a  fixed  k  and  j  a  k+1,  k+2,....  where  the  initial  condition  is  x(k|k). 

The  gain  matrix  B  is  given  by 

j-1 

B(j)  =  x  A(i)  (38) 

i=k 

and 

A(i)  »  P(ij  i)0t(i+l;i)P“1(i+lj  i)  (39) 

The  fixed  point  smoothing  error  covariance  matrix  is  given  by 

P(k|j)  =  P(k|  j-1)  +  B(j)(P(j|j)-P(j|j-l))Bt(j)  (40) 

An  alternate  formulation  of  this  equation  is  possible  but  was  not  used  in 
this  study  (Ref  14:231). 

Reverse  Integration 

Once  the  final  filtered  estimate  has  been  obtained  using  the  equations 
(30)  to  (33) ,  the  final  state  can  be  reverse  integrated  using  the  non¬ 
linear  state  equations  (13)  so  as  to  obtain  the  state  at  any  previous  time 
as  reflected  by  the  state  at  the  final  measurement.  A  point  of  special 
interest  is  the  initial  point  since  reverse  integration  to  this  point 
allows  comparison  with  the  results  of  the  Weighted  Least  Squares  filter 
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presented  in  chapter  IV. 

€ 

Class  I  Filter  Modification 

Class  I  Filter  Modification  was  employed  throughout  this  study.  This 
subject  was  first  covered  in  Reference  19.  This  attempt  to  account  to 
some  extent  for  computer  truncation  and  round-off  error  as  well  as  system 
non-linearities  not  adequately  represented  by  the  system  model  takes  the 
form: 

P(k  +  l|k)  =  0(k  +  l,k)P(k|k)0T(k  +  l,k)  + 

q 


The  value  of  q  is  a  function  of  the  amount  of  non-linearity  and  round-off 
anticipated.  For  most  of  the  study  a  value  for 

q  o  1  x  10-12  (42) 

was  used  although  for  purposes  of  analyses  it  was  varied  at  certain  points. 
This  will  be  further  discussed  in  Chapter  V  when  the  simulated  orbit 
analysis  is  performed. 


t 
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Weighted  Least  Square  Error  Filterinc 


Statistical  filtering  or  least  square  error  filtering  assumes  a  cost 
function  (Ref  18:16)  indicative  of  the  desired  outcome  of  the  problem. 

The  nominal  trajectory  is  to  be  chosen  such  that  when  it  is  compared  to 
the  values  of  the  measured  trajectory  the  difference  of  the  squared  error 
is  minimized.  If  the  errors  in  certain  parameters  are  weighted  more  than 
the  errors  in  other  parameters  than  the  method  is  called  Weighted  Least 
Squares  Filtering.  The  statistical  filter  then  provides  both  a  smoothed 
solution  and  a  filter  solution  simultaneously  when  the  problem  is  solved 
since  the  method  is  necessarily  interested  in  the  errors  found  along  the 
entire  trajectory. 

Unlike  the  probabilistic  or  Kalman  filter,  the  Least  Squares  Filter 
is  a  batch  filter  which  requires  the  simultaneous  application  of  all  data 
to  be  processed.  The  addition  of  one  new  data  point  under  the  batch 
process  will  require  the  rerunning  of  the  entire  "batch"  of  data.  This  is 
the  chief  disadvantage  of  the  statistical  filter;  however,  methods  are 
available  to  transform  this  batch  processor  into  a  recursive  processor 
through  changes  in  the  cost  function. 

A  cost  function  which  expresses  the  least  square  error  concept  stated 
above  is  given  by: 

j  =  |  lx  -  y(x0)||J  3  (i-yQ^Acx-yUo))  (43) 

This  quadratic  cost  function  must  be  convex  and  therefore  has  a  minimum. 

Using  minimization  techniques  common  to  optimal  control  and  the 
calculus  of  variations  it  can  be  shown  that 

AXq  =  (ATWA)‘1ATW  Ay  (44) 

will  cause  this  cost  function  to  be  minimized,  where 
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A  =  the  measurement  sensitivity  matrix 

=  the  required  change  in  the  nominal  trajectory 
and 

W  =  the  weighting  matrix 

The  Measurement  sensitivity  matrix.  A,  can  be  determined  by  use  of 
the  measurement  matrix  M  previously  discussed  in  Chapter  IV  and  Appendix  D. 

v<w ' 
wv 

M20(t2,to)  (45) 

V'VV 

where  0  is  the  state  transition  matrix  found  from  the  solution  of  the 
linear  matrix  differential  equation 

0(t,to)  -  F  0(t,to)  (46) 

Unfortunately  due  to  round-off  and  truncation  error  as  well  as 
modeling  non-linearities  the  value  found  for  AXg  will  not  cause  the  cost 
function,  J,  to  reach  its  minimum  without  repeated  iteration  of  this  pro¬ 
cess.  Therefore,  some  criteria  must  be  set  for  determining  convergence 
to  the  "true"  state  xQ. 

The  criteria  chosen  involve  the  residuals.  They  are  forced  to  be 
less  than  some  arbitrarily  small  value.  Care  must  be  exercised  to  make 
this  value  realizable  since  if  it  is  made  too  small,  then  many  iterations 
will  be  required  to  obtain  convergence.  If,  on  the  other  hand,  it  is  not 
small  enough  then  the  state  Xq  will  not  be  accurately  determined. 

For  this  study  a  value  for  the  epsilon  of  convergence  equal  to 
2  x  10”^  was  found  to  be  suitable. 
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If  we  define 


then  if  either 


or 


I  UyP|  |  «  ||Ay  -  AAX  |1 

Wo  W 

r  _  l|Ay||w  -  ! I AyP| |W  _  _ 

^1  . .  .  >  c. 

I  Ml, 

c2  *  I  l*y|  lw  -  |  |iyP|  |  <  c. 


were  true  then  the  data  was  said  to  have  converged. 


(47) 


(48) 


(49) 
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To  verify  filter  operations  and  to  provide  a  basis  for  analysis  of 
actual  radar  data,  a  simulated  orbit  analysis  was  first  performed  using 
both  the  WLS  and  Kalman  Filtering  methods  on  an  orbit  having  known  orbital 
parameters. 


The  Simulated  Orbit 

The  following  parameters  were  chosen  for  the  simulated  orbit: 
tQ:  Day  298  Hour  02  Minute  39  Second  19.466 

6  :  Day  298  Hour  02  Minute  12  Second  18.503  (Ref  Appendix  F) 

go 

Serai-Major  Axis  :  1.161913  earth  radii  =  aQ 

Eccentricity  :  .06539071  =  eQ 

Angle  of  Inclination:  54.2538°  =  i 
Line  of  Nodes  :  199.6287°  = 

Argument  of  Perigee  :  29.9570°  a  u)Q 

Period  :  105.8204  minutes  =  PQ 

Height  at  Perigee  :  548.174  Kilometers 

The  period  and  height  at  perigee  are  functions  of  the  other  orbital 
elements  and  are  presented  only  since  they  are  often  parameters  of  inter¬ 
est. 

Care  must  be  exercised  to  associate  the  time,  tQ,  with  these  orbital 
parameters  since  the  perturbation  effects  caused  by  the  oblate  earth,  non- 
uniform  gravitational  field,  and  the  time-varying  density  models  will 
cause  these  parameters  to  vary  as  the  satellite  orbits  the  earth. 


Simulated  Orbit  Data  Generation.  With  the  above  orbital  parameters 
as  the  initialization  point,  a  two  pass  orbit  was  generated  making  use 
of  a  single  simulated  radar  whose  parameters  are  given  in  Table  I. 
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Table  I 

Theoretical  Radar  Parameters 

Position 

Latitude  =  52.73267°N 

Longitude  =  174.1023CE 
Height  =  93.  Meters 

Measurement  Biases 
Range  =0.0 
Range  Rate  =0.0 
Azimuth  =  0.0 
Elevation  =  0.0 

Measurement  Sigmas 
Range  =  1000  meters 
Range =1  meter/second 
Azimuth  =  .02  degrees 
Elevation  =  .02  degrees 


Gaussian  Distribution  Data  Analysis 
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The  height:  of  the  radar  is  measured  with  respect  to  the  reference  ellip¬ 
soid. 

In  order  to  obtain  a  noisy  observation,  true  measurements  of  azimuth, 
elevation,  range,  and  range  rate  were  contaminated  with  computer  gener¬ 
ated  noise.  By  the  repeated  addition  of  points  selected  from  uniform  dis¬ 
tributions  an  approximately  normal  distribution  can  be  obtained  from 

12 

N(0,1)  s  £  x.  -6  (50) 

i=l  1 

where  the  x^  are  random];*  selected  from  a  uniform  distribution  between  0 
and  1.  This  distribution  is  known  to  be  a  truncated  distribution  having 
at  best  a  variance  of  six  sigma.  However,  for  the  orbital  problem  at  hand 
variations  outside  this  range  are  unlikely  and  the  distribution  is  there¬ 
fore  accepted  as  adequate. 

A  study  of  the  generated  noise  distributions  was  made  using  the  shape 
of  the  generated  distribution  curve  as  a  rough  measure  of  the  density  of 
the  distribution.  Figure  3  shows  a  plot  of  a  distribution  composed  of  500 
sample  points.  Other  plots  of  various  numbers  of  points  are  presented  in 
Appendix  E.  Other  basic  statistical  descriptors  of  the  samples  are  shown 
in  Table  II.  The  non- zero  mean  shown  in  Table  II  presents  a  small  problem 
to  the  filters  since  the  versions  of  both  which  are  programmed  for  this 
study  expect  zero  mean  noises  as  would  be  anticipated  in  the  actual  orbital 
determination  problem.  Small  biasing  effects  should  therefore  be  antici¬ 
pated  in  the  results.  The  samples  indicated  were  also  subjected  to  Kolo- 
Eogorov -Smirnov  Test  with  a  significance  level  of  .20,  and  all  of  the  sam¬ 
ples  were  accepted  as  possibly  being  from  a  N(0,1)  distribution. 

On  the  basis  of  the  distributions  studied,  it  was  determined  that  a 
single  run  consisting  of  two  sets  of  fifty-five  measurements  each  of  azi¬ 
muth,  elevation,  range,  and  range  rate  could  not  present  a  truly  Gaussian 
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noise  to  the  filters.  Therefore,  a  series  of  25  runs  against  the  same 
theoretical  parameters  but  contaminated  with  different  sequences  of  noise 
were  studied  with  the  statistical  average  of  the  parameters  determined 
being  accepted  as  the  final  determination  of  the  filters  in  the  theoreti¬ 
cal  study. 

Orbital  Perturbations.  Figures  2  through  6  show  the  true  orbital 
parameters  as  determined  using  the  previously  described  models  of  the 
oblate  earth,  non-uniform  gravitational  field,  and  time-varying  atmosphere. 
The  large  tick  on  the  abscissa  indicates  the  passage  of  one  revolution. 

Table  III  summarizes  the  effects  noted  on  these  graphs.  Generally,  these 
perturbations  can  be  divided  into  two  types,  secular  and  those  due  to  the 
oblate  earth  effects  (Ref  21:208-241).  Appendix  G  notes  how  the  state  vector 
may  be  transformed  into  the  classical  orbital  elements  graphed. 
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Figure  4.  True  Simulated  Semi-Major  Axis 


Figure  7.  True  Simulated  Line  of  Nodes 
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Table  III 


Orbital  Perturbation  Analysis 


Parameter 

Max 

Min 

Cycles 

Secular 

Change 

(Deg/Rev) 

Semi -Major 

Axis  (E.R. ) 

1.1637 

1.1616 

2 

none 

Eccentricity 

.0667 

.0654 

3 

none 

Period  (Min) 

105.78 

106.07 

2 

none 

Angle  of 

Inclination 

(Deg) 

54.282 

54.247 

2 

none 

Line  of  Nodes 
(Deg) 

199.65 

199.40 

2 

-.25° 

Arg.  of 

Perigee  (Deg) 

30.435 

29.443 

3 

+  .15° 
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The  graphs  and  analytic  results  from  the  integration  of  the  state  equa¬ 
tions  yield  the  following  compatible  results: 

LSI  =  199.347°  -  199.602°  =  -  .255  deg/rev  (55) 

Aw  =  +  29.700°  -  29.551°  =  .149  deg/rev  (56) 


Oblate  Earth  Perturbations.  Since  this  study  is  confined  to  low  eccen¬ 
tricity  orbits  ranging  from  approximately  circular  (0)  to  a  value  of  about 
.124,  certain  approximations  to  the  general  equations  given  in  Reference 
21,  page  220,  can  be  made.  Basically  all  terms  which  contain  the  eccen¬ 
tricity  at  t  and  powers  of  it  are  deleted  as  being  small  with  respect  to 
1.  Thus  these  equations  can  be  reduced  to: 


a  s  K_  +  - if|e -  sin2i_cos2u„ 


(57) 


2a0a-e02)3 


e  2  k  ♦  -..y^e 

C 


2a02(l-e02)2 


^ (1-3/2  sin2i0)cosf°  +  i  sin2i0cos(2wo+f°) 
♦  7/12  sin2iocos(2wo+3f0)  ) 


3J9ao2sin2i_ 
i  s  k.  +  — l-S. - ° 

Ha  2 n  o  2>2 

8ao  U-eo  ) 


(  cos2u0) 


,  .  K„  ■  -  (2f°  -  si»2u0) 


4a02(l-e02)2 


“  E  ’  4»^(l-e0V  (  (4'5sin2io>f°  *  ~  (2-3sin2i„)sin  fC 


(59) 

(60) 


(61) 


-  ^  sin2iQsin(2w0  +  f°)  +  sin2iosin(2w0  +  3f0)) 
where  secular  perturbation;  terms  have  been  included  in  fi  and  w  and  where 


uQ  =  wQ  ♦  f 


and 


.  /  H 

dt  \l  aQ3(l-e2)3 
is  the  true  anomaly  at  tQ. 


(1  +  eQccs 


and  f® 


(62) 
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Since  the  true  anamaly,  f,  must  progress  through  360  physical  degrees 
per  revolution  or  one  cycle,  the  number  of  cycles  in  the  predominant  fre¬ 
quency  component  of  the  perturbations  can  be  predicted  from  the  preceeding 
equations.  Thus  it  can  be  anticipated  that  the  semi-major  axis,  period, 
angle  of  inclination,  height  at  perigee,  and  the  line  of  nodes  will  exper¬ 
ience  two  cycles  per  revolution  while  the  eccentricity  and  the  argument  of 
perigee  will  experience  three  cycles.  These  conclusions  are  borne  out  by 
the  figures  in  this  chapter.  It  should  also  be  anticipated  that  the  argu¬ 
ment  of  perigee  will  become  extremely  difficult  to  determine  as  the  eccen¬ 
tricity  of  the  orbit  approaches  zero.  If  the  eccentricity  and  argument 
of  perigee  become  extremely  difficult  to  calculate,  a  change  in  orbital 
parameters  can  be  made  as  detailed  in  Reference  21,  page  194.  However, 
in  this  study  no  change  in  parameters  will  be  made  since  the  Ent  analyses 
are  based  on  parameters  of  the  type  determined  in  Appendix  G. 

Filter  Operation.  Both  the  Kalman  Filter  and  WLS  filter  presented 
highly  accurate  estimates  of  the  orbital  parameters.  Comparison  of  the 
estimated  parameters  occurs  in  Tables  IV  and  V.  The  errors  in  the  out¬ 
puts  of  the  Kalman  filter  as  a  function  of  the  number  of  measurements  are 
presented  in  Figures  7  through  20.  It  must  be  noted  that  the  points 
plotted  in  these  figures  are  at  observations  1,2,3,4,5,10,  and  every  fifth 
observation  thereafter  during  the  first  pass  and  at  observations  56,57,58, 
59,60,65  and  every  fifth  observation  thereafter  in  the  second  pass.  Thus 
general  trends  are  being  displayed  with  a  change  of  scale  usually  occurring 
on  the  second  graph  in  each  set  so  that  the  data  may  be  more  accurately 
presented.  Generally,  the  data  converges  toward  a  zero-mean  error.  How¬ 
ever,  a  bias  of  .02°  is  seen  in  the  line  of  nodes  parameter.  Previous 
discussion  of  the  noise  distribution  involved  may  account  for  this  bias. 
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Non-linearity  due  to  initial  condition  errors  in  the  early  phases  of  the 
filtering  processes  also  undoubtedly  contributed  to  all  of  the  residual 
errors.  Both  the  values  shown  and  the  errors  graphed  are  the  results 
of  the  25  run  statistical  averages  discussed  previously.  All  parameters  - 
determined  are  considered  sufficiently  accurate  to  enable  future  position 
prediction.  The  values  determined  for  the  orbital  elements  by  the  Kalman 
filter  are  shown  in  figures  21  through  25. 
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Table  IV 


Kalman  Filter  Analysis* 


Number  of 
Observations 

20 

30 

40 

55 

75 

110 

Error  In 

Semi-Major 
Axis  (E.R.) 

-.000078 

.000014 

.000008 

-.000021 

.000073 

.000051 

Eccentricity 

-.0000463 

.0000205 

.0000156 

.0000130 

.0000597 

.0000287 

Angle  of 

Inclination 

(Deg) 

-.0004 

-.0003 

-.0002 

-.0002 

-.0005 

.0001 

Line  of 

Nodes  (Deg) 

.232 

.0246 

.0246 

.0248 

.0347 

.0246 

Angle  at 
Perigee  (Deg) 

-.0380 

-.0171 

-.0169 

-.0128 

.0048 

.0188 

Height  at 
Perigee  (Km) 

-.124 

-.072 

-.070 

-.044 

-.007 

-.069 

Period 

-.010 

.002 

.001 

.001 

.010 

.007 

•This  data  is  also  presented  graphically  in  figures  7  through  20 
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Table  V 

WLS  Filter  Analysis 
Epoch  time  =  tQ  =  0.0  seconds 


Number  of 
Observations 

20 

30 

40 

55 

75* 

Error  In 

Semi-Major 

Axis  (E.R.) 

.000014 

.000046 

.000018 

.000009 

.000049 

Eccentricity 

.0000146 

.0000438 

.0000260 

.0000180 

.0000836 

Angle  of 

Inclination  (Deg) 

.0000 

.0002 

.0002 

.0001 

.0044 

Line  of 

Nodes  (Deg) 

-.0047 

-.0038 

-.0040 

-.0041 

-.0015 

Angle  at 

Perigee  (Deg) 

.0003 

-.0038 

-.0125 

-.0139 

-.0570 

Height  at 

Perigee  (Km) 

1 

• 

O 

K> 

-.050 

-.084 

-.078 

.323 

Period  (Min) 

.0019 

.0063 

.0025 

.0012 

.0067 

Iterations 

3 

3 

3 

3 

7 

*The  75  point  run 

is  a  single 

run  and  not 

the  average 

of  25  runs 

as  are 

the  others.  This  was  necessitated  by  the  large  amount  of  machine  time 
required  by  repeated  integrating  the  state  equations  around  the  orbit  with 
different  starting  conditions.  (CP  Time,  Single  Run  75  Points  =  484  sec). 
However,  convergence  on  pass  two  was  achieved  after  7  iterations. 


flLMAN  ERRORS 


KALMAN  ERRORS  PASS 


Figure  10,  Semi -Major  Axis  Error  Pass 
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Figure  11.  Eccentricity  Error  Pass 
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Figure  13.  Angle  of  Inclination  Error  Pass 
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Figure  14.  Angle  of  Inclination  Error  Pass 
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Kalman  Determined  Semi -Major  Axis,  Simulated  Orbit 


THEORETICAL  DATA 


Figure  24.  Kalman  Determined  Eccentricity,  Simulated  Orbit 


TICAL  DATA 


Figure  25.  Kalman  Determined  Angle  of  Inclination,  Simulated  Orbit 
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Figure  27.  Kalman  Determined  Period,  Simulated  Orbit 
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Both  types  of  smoothing  discussed  in  Chapter  III 


were  applied  to  the  theoretical  data. 


Fixed  Interval  Smoothing.  Fixed  interval  smoothing  failed  when 
applied  to  the  simulated  orbit  data.  After  processing  all  data  points  in 
the  first  pass  with  the  filtering  algorithm,  the  smoother  was  applied  to 
attempt  to  smooth  back  to  point  SIS.  However,  at  point  #47  the  smoothing 
covariance  matrix  lost  symmetry  and  soon  thereafter,  negative  terms  began 


to  appear  on  the  diagonal.  Because  of  these  effects,  this  method  of 


smoothing  was  abandoned. 


Fixed  Point  Smoothing.  The  fixed  point  smoother  performed  extremely 
well  when  applied  during  the  first  pass  of  the  theoretical  data.  Table  VI 
shows  some  comparative  results.  As  can  be  seen  some  of  the  best  and  some 
of  the  worst  errors  occur  with  this  smoother.  For  those  terms  which  do  not 
experience  secular  variations,  the  smoother  works  with  extreme  accuracy. 
However,  for  the  line  of  nodes  and  the  argument  of  perigee  large  errors 
occur. 

Less  successful  results  were  obtained  from  this  smoother  during  second 
pass  applications.  The  probable  cause  for  this  change  in  performance  was 
the  change  in  the  behavior  of  the  filtering  covariance  matrix  in  the  second 
pass.  This  point  will  be  discussed  thoroughly  in  the  following  section 
concerning  the  covariances  encountered  using  the  simulated  orbit  data.  No 
problems  were  encountered  with  loss  of  symmetry  or  negative  diagonal  terms 
while  using  the  fixed  point  smoother. 

All  smoothing  performed  during  the  simulated  orbit  study  was  done 
at  point  #15  or  with  respect  to  it.  The  choice  of  this  point  was  arbitrary 
although  it  was  chosen  with  the  intention  of  allowing  the  filtering  algorithm 
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Table  VI 

Kalman  Smoother  Analysis 


Point  #  15  is  being  smoothed.  t=84  seconds. 


Number  of 
Observations 

20 

30 

40 

55 

Errors  In 

Semi-Major 

Axis  (E.R. ) 

-.000086 

.000005 

.000001 

.000000 

Eccentricity 

.0000479 

-.0000196 

.0000149 

.0000121 

Angle  of 

Inclination  (Deg) 

.0007 

,0004 

.0003 

,0002 

Line  of 

Nodes  (Deg) 

.1229 

.3749 

.6257 

1.0020 

Angle  at 

Perigee  (Deg) 

.0462 

.0248 

.0245 

.0203 

Height  at 

Perigee  (Km) 

.193 

.113 

.063 

.090 

Period  (Min) 

.012 

.001 

.000 

.000 

C 
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to  have  a  chance  to  settle  towards  a  steady  state  before  the  smoother  was 
applied. 

Covariance  Matrix  Behavior.  The  behavior  of  the  covariance  matrix 
during  this  study  was  typical  of  those  matrices  described  in  theory  (Ref 
15) .  Starting  with  an  intentionally  large  and  diagonal  covariance  at  t 
the  covariance  matrix  continues  a  rapid  constant  descent  until  a  "steady 
state"  is  achieved.  Figures  26  through  29  show  the  behavior  of  the  Pxx 
or  Pjj  element,  the  P^j  or  P44  element,  as  well  as  the  Position  and  Veloc¬ 
ity  Standard  Deviation  which  are  the  square  root  of  the  sum  of  the  squares 
of  the  appropriate  covariance  elements: 

o  o  2  1/2 

Position  Standard  Deviation  =  (Pxx  +  Pyy  +  Pzz  )  '  (63) 

2  2  21/2 

Velocity  Standard  Deviation  =  (P**  +  P;  *  +  P;;  )  (64) 

AX  yy 

Initial  Covariance  Matrix,  P(0).  The  initial  covariance  matrix  was 
chosen  after  experimentation  with  various  values.  Values  as  large  as 
1G30  for  position  variance  and  10*5  for  velocity  covariance  were  tried. 

These  values  converged  rapidly  to  the  same  steady  state  values  as  those 
obtained  with  the  values  actually  used,  1010  and  106.  Values  between  these 

limits  were  also  tried  with  convergence  to  the  same  values  occurring. 

5  3 

However,  it  was  determined  that  if  small  values  were  chosen,  10  and  10  , 
for  example,  then  convergence  to  the  same  values  did  not  occur  and  in  fact 
a  steady  state  was  never  achieved  during  the  allowed  observation  period  of 
110  measurements.  Additionally,  physical  significance  can  be  attributed 
to  values  of  10*®  for  position  and  10®  for  velocity.  These  values  represent 
one  sigma  values  of  10  and  10*  for  position  and  velocity  respectively. 
Because  of  the  manner  in  which  the  initial  state  x(0)  is  chosen,  little 
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can  truly  be  said  about  the  position  and  velocity  of  the  vehicle  other 
than  that  orbital  speed  was  attained.  Since  the  position  and  velocity 
are  of  the  order  of  magnitude  of  10  and  10  meters  and  meters/second 
respectively,  this  choice  of  sigma  implies  that  very  little  statistical 
confidence  should  be  placed  in  the  initial  values. 

Experimentation  with  the  Kalman  Gain,  K.  During  the  study  the 
Kalman  matrix  was  varied  by  premultiplication  by  a  constant  and  by  pre¬ 
multiplication  of  columns  of  elements  by  different  constants.  Theory 
shows  that  the  Kalman  matrix  as  derived  by  Kalman  is  the  optimal  which 
can  be  applied  to  a  set  of  data  to  obtain  estimates  of  the  states.  Ex¬ 
perimentation  confirms  this  result  and,  in  fact,  indicated  that  the  Kalman 
matrix  is  extremely  sensitive  to  manipulation. 

Multiplication  by  constants  varying  from  .8  to  3.  was  attempted. 
Multiplication  by  a  constant  less  than  one  provided  slower  convergence  to 
the  values  achieved  than  when  the  unaltered  K  matrix  is  used.  Multipli¬ 
cation  by  a  constant  greater  than  one  tended  to  induce  oscilations  in  the 
values  of  the  elements  until  when  the  constant  was  increased  greater  than 
2,  filter  divergence  occured. 

Premultiplication  by  different  patterns  of  constants  produced  no 
acceptable  effects  for  the  orbital  determination  problem. 

Experimentation  with  the  value  of  P.  The  value  of  the  covariance 
matrix  P  was  varied  by  the  addition  of  a  constant  term  to  the  values  for 
the  velocity  covariances  at  every  measurement  update  in  accordance  with 
recommendations  for  compensation  for  non-linearity  and  covariances  which 
approach  zero  (Ref  10:305-307).  This  effective  noise  addition  produced 
no  desirable  results  over  the  number  of  observations  used.  Oscillations 
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about  the  value  obtained  using  the  unaltered  covariance  matrix  tended  to 
occur „ 

Experimentation  with  the  value  of  Q.  The  value  of  Q  chosen  for  use 
with  Level  One  Modification  (Ref  19)  was  1.  x  10  .  The  number  of  sig¬ 

nificant  digits  carried  by  the  CDC  6600  in  single  precision  arithmetic  is 

14.  Thus  to  be  exact  theory  requires  a  value  of  1.  x  10  .  The  lower 

-12 

value  of  10  assists  in  accounting  for  additional  non-linearities  in 

-20 

the  system.  However,  it  should  be  noted  that  changes  to  values  of  10 
and  zero  were  made  to  the  modification  with  little  noticeable  effect.  The 
large  capacity  and  capability  of  the  CDC  6600  probably  precludes  necessity 
for  this  type  of  modification.  The  between  pass  covariance  matrix  reached 
a  maximum  of  105  for  the  diagonal  position  terms  and  103  for  the  diagonal 
velocity  terms  prior  to  update  using  the  first  measurement  of  the  second 
pass.  No  tendency  toward  divergence  was  observed. 

Pointing  Information  from  Prediction 

The  use  of  prediction  capabilities  of  the  filters  is  of  key  importance 
in  the  pointing  of  downrange  radars  in  order  to  obtain  additional  tracking 
data  for  orbital  refinement  and  also  to  enable  possible  rendezvous  or 
interception.  Table  VII  shows  the  capability  of  the  Kalman  filter  to  pre¬ 
dict  future  measurements.  Two  items  are  important  when  comparing  these 
numbers.  First  the  time  of  observation  and  prediction  may  vary  slightly 
and  second  the  measurement  is,  of  course,  corrupted  by  noise  while  the  fil¬ 
ter  parameters  are  best  estimates. 
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Table  VII 

Theoretical  Interpass  Pointing  jta  Kalman  Filter 
Single  Run  (randomly  chosen) 


Pass 

Range 

(meters) 

Range  Rate 
(meters/second) 

Azimuth 

(Deg) 

Elevation 

(Deg) 

1-2  Calculated* 

1986441 

-6530 

259,437 

10.327 

1-2  Measured 

1982343 

-6526 

259.437 

10.384 

^Calculations  are  performed  using  the  equations  found  in  Appendix  C. 


( 
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Vic  Actual  Data:  Ent  Set  1 

Data.  Data  for  approximately  three  days  of  the  orbit  of  satellite 
6096  was  supplied  for  this  study  by  the  1st  Aerospace  CON  Squadron,  Ent 
AFB,  Colorado.  This  data  was  processed  at  Ent  using  a  special  perturba- 
tional  analysis  and  the  results  were  provided  for  comparison  in  this 
study.  The  Ent  analysis  used  approximately  660  observations  of  the  ve¬ 
hicle,  rejecting  about  40  observations  for  failure  to  fit  the  predicted 
statistics  of  the  problem.  The  spiral  decay  differential  correction  pro¬ 
gram  converged  providing  orbital  elements  at  an  epoch  time  corresponding 
to  the  final  measurement  on  the  third  day. 

From  the  supplied  data,  115  observations  were  selected  to  be  pro¬ 
cessed  by  the  Kalman  Filtering  program.  These  observations  were  chosen 
because  they  represent  a  three  orbit  two  radar  problem.  Nothing  was 
known  of  the  characteristics  of  any  particular  observation  at  the  time  of 
selection.  The  Kalman  program  produced  orbital  elements  at  the  time  of 
its  final  measurement  which  occurred  early  in  the  first  day  of  orbit.  To 
obtain  values  of  the  orbital  elements  at  other  times  a  reverse  integration 
to  the  beginning  was  accomplished.  To  predict  values,  the  fact  that  the 
elements  remain  periodic  and  have  known  secular  variations  was  used.  This 
allows  comparison  of  the  determined  elements  with  elements  generated  by 
the  other  programs  at  any  time. 

The  first  80  of  the  same  115  data  points  were  used  to  run  the  WLS  fil¬ 
ter.  A  forward  integration  was  performed  to  project  the  determined  ele¬ 
ments  to  other  times  for  comparison  with  ENT  and  Kalman  filtering  data. 
Table  VIII  shows  the  specifics  of  the  data  used. 
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Table  VIII 
Orbit  Data 

Satellite  #  6096  Day  190  Year  72 

0  =  19  hr  4  min  23.4378  sec  (Ref  Appendix  F) 

Ko 

Start  GMT  =  20  hr  41  min  8.839  sec 


Pass 

Observation 

Interval 

Number  of 
Observations 

Duration 

Pass 

Min  Sec 

Station 

Duratic 
of  Gap 
Min  Sec 

la 

10  sec 

22 

3:30 

359 

lb* 

10  sec 

25 

4:10 

359 

1:15.7 

2 

10  sec 

33 

5:30 

359 

85:29.5 

3 

6  sec 

35 

3:30 

345 

91:30.3 

*Pass  lb  is  a  continuation  of  pass  la  which  occurred  due  to  a  radar  dropout 
of  1  minute/15  seconds. 


( 
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Radar  Characteristics.  The  parameters  or  characteristics  of  the 
radars  used  are  shown  in  Table  IX.  The  similarities  between  these  proper¬ 
ties  and  those  of  the  simulated  radar  discussed  in  Chapter  V  are  noted. 
Passes  la,  lb,  and  2  are  obtained  from  radar  359  while  pass  3  comes  from 
radar  345. 

Comparison  of  Orbital  Parameters.  The  results  of  the  Kalman  filter 
and  the  WLS  filter  projected  forward  to  the  epoch  time  of  the  ENT  data 
are  compared  in  Table  X.  From  this  table  a  high  degree  of  compatibility 
is  noted  in  all  parameters  except  the  argument  of  perigee.  Since  this 
orbit  is  extremely  circular  in  nature,  the  very  term  perigee  tends  to  lose 
meaning  and  thus  so  does  the  parameter  argument  of  perigee.  As  a  result 
this  parameter  is  seen  to  vary  from  290°  to  355°.  The  eccentricity  appears 
to  have  a  larger  error  than  the  other  parameters  but  again  it  is  noted  that 
an  eccentricity  such  as  those  obtained  merely  indicates  that  the  orbit  is 
circular.  A  maximum  difference  of  .007  minutes  is  noted  between  the  WLS 
period  and  the  ENT  period.  This  comes  to  .42  seconds  which  is  certainly 
an  acceptable  error  considering  the  small  amount  of  data  processed.  Graphs 
of  the  Kalman  determined  orbital  elements  are  found  in  figures  30  through 
34.  These  graphs  exhibit  the  same  characteristics  noted  in  the  theoreti¬ 
cal  study  although  the  parameters  have  assumed  a  more  symmetric  shape  due 
to  the  polar  circular  nature  of  the  orbit. 


73 


-r. 


++*&»**  i* 


GGC/EE/73-13 


Table  IX 

Radar  Parameters  Ent  Set  1 
Position 


Sensor  8  359 

345 

Latitude 

64.291  N 

52.73267  N 

Longitude 

149.19318  ft 

185.8976  E 

Height 

199. 

86.  meters 

Sensor  Biases 

Range 

105.9 

218.3  meters 

Range  Rate 

1.0 

1.0  meters/second 

Azimuth 

-.037 

-.067  degrees 

Elevation 

.023 

.054  degrees 

Sensor  Sigmas 

Range 

819. 

929.  meters 

Range  Rate 

24. 

18.  meters/second 

Azimuth 

.474 

.259  degrees 

Elevation 

.331 

.268  degrees 

1 
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Table  X 


Ent  Data 

Set  I 

Results 

Epoch  Time: 

Day:  192 

Hour: 

20  Min:  01 

Sec:  23.325 

Kalman 

WLS 

Element 

ENT 

Predicted 

Predicted 

Semi -Major 

Axis 

1.076376 

1.076902 

1.076808 

Eccentricity 

.000915 

.00160 

,00141 

Angle  of 
Inclination 

96.149° 

96.152° 

96.154° 

Line  of  Nodes 

N/A* 

248.367° 

248.369° 

Arg.  of  Perigee 

290.2294° 

354.889° 

345.234° 

Period 

94.415  min 

94.420  min 

<0 

• 

O 

00 

y 

u. 

Height  at 

Perigee 

480  km 

480  km 

480  km 

#  of  Observations 
Used  to  Obtain 
Results 

660 

115 

80 

Radius  (m) 

6871522 

6873301 

6873552 

Velocity  (m/sec) 

7612.8 

7612.7 

7612.1 

*Not  Available 
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Figure  34.  Angle  of  Inclination  Ent  Set 
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Kalman  Smoothing.  The  smoothing  capability  exhibited  by  the  fixed 
point  smoother  using  simulated  orbit  data  was  not  as  effective  when  applied 
to  true  data.  Primary  causes  for  this  are  the  loss  of  a  white  Gaussian 
distribution  due  to  a  smaller  number  of  observations  applied  to  the  smooth¬ 
ing  algorithm  as  well  as  the  inability  of  the  filter  to  stabilize  the 
orbital  parameters  as  quickly  using  the  real  data.  However,  stability 
was  achieved  during  passes  2  &  3  with  the  results  comparing  favorably 
with  the  Kalman  reverse  integration  as  shown  in  Table  XI.  In  order  to 
obtain  a  maximum  number  of  measurements  the  point  being  smoothed  was  the 
first  observation  point  for  each  set  of  data.  This  was  done  at  the  risk 
and  with  the  understanding  that  the  desired  stability  obtained  by  waiting 
till  the  fifteenth  point  as  in  the  theoretical  study  would  be  lost.  Al¬ 
though  satisfactory  results  were  not  obtained  for  first  pass  data,  no 
extremely  adverse  effects  were  noted  and  the  smoother  was  gradually 
bringing  the  unstable  first  observation  parameters  toward  the  actual 
parameters  although  hampered  severely  by  a  lack  of  measurements  in  each 
of  the  first  two  data  sets  (la  and  lb).  The  secular  varying  terras,  line 
of  nodes  and  argument  of  perigee,  are  again  noted  as  failing  in  all  cases 
to  reach  an  acceptable  value.  Smoothing  at  points  five  and  fifteen  was 
also  performed  with  no  appreciable  difference  in  the  results. 
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Table  XI 


Smoothing  Data  at  First  Observations 


Ent  Set  1 


Kalman 

Reversed 

WLS 

Smoother 

Pass  la 

t=0. ,  22  observations  used  for  smoothing. 

Semi-Maj or 

Axis  (E.R.) 

1.07638 

1.07629 

1. 07342 

Eccentricity 

.00128 

.00126 

.00374 

Angle  of 

Inclination  (Deg) 

96.153 

96.155 

96.160 

Line  of  Nodes  (Deg) 

246.739 

246.745 

247.619 

Arg.  of  Perigee 
(Deg) 

332.050 

246.745 

288.591 

Period  (Min) 

94.352 

94.340 

93.963 

Pass  lb 

t=285.7.  25  observations  used  for  smoothing. 

Semi -Major 

Axis  (E.R.) 

1.07695 

1.07676 

1.07701 

Eccentricity 

.00162 

.00141 

.00171 

Angle  of 

Inclination  (Deg) 

96.151 

96.153 

96.162 

Line  of 

Nodes  (Deg) 

246.746 

246.749 

247.733 

Arg.  of  Perigee 
(Deg) 

3.711 

347.809 

6.378 

Period  (Min) 

94.429 

94.395 

94.541 
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Smoothing  Data  at  First  Observations  (Cont.) 


Pass  2 

t=5655.2.  33  observations  used  for  smoothing.* 


Semi -Major 

Axis  (E.R. ) 

1.07636 

1.07627 

1.07639 

Eccentricity 

.00126 

.00119 

.00116 

Angle  of 

Inclination  (Deg) 

96.153 

96.156 

96.154 

Line  of  Nodes  (Deg) 

246.795 

246.799 

248.136 

Arg.  of 
t‘erigee  (Qeg) 

336.732 

315.589 

313.768 

Period  (Min) 

94.348 

94.336 

94.353 

Pass  3 

t=11465.5.  35  observations  used 

for  smoothing.* 

Semi -Major 

Axis  (E.R.) 

1.07655 

1.07645 

1.07655 

Eccentricity 

.00149 

.00139 

.00150 

Angle  of 

Inclination  (Deg) 

96.153 

96.155 

96.153 

Line  of  Nodes  (Deg) 

246.851 

246.855 

247.704 

Arg.  of 

Perigee  (Deg) 

344.000 

332.569 

344.984 

Period  (Min) 

94.374 

94.361 

94,371 

*A  steady  state  was  achieved. 
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Interpass  Pointing  Data.  Table  XII  shows  the  values  for  the  measure¬ 
ments  and  the  predicted  values  of  the  range,  range  rate,  azimuth  and  ele¬ 
vation.  These  results  are  extremely  similar  to  those  obtained  in  the 
theoretical  case  and  certainly  are  adequate  to  provide  pointing  informa¬ 
tion. 


Filter  Covariance.  Graphs  of  the  position  and  velocity  covariance 
of  the  x  parameter  are  shown  in  figures  35  and  36.  These  graphs  exhibit 
only  slightly  differing  characteristics  from  those  found  in  the  theoreti¬ 
cal  case.  The  time  lag  between  pass  la  and  lb  is  seen  to  be  small  enough 
that  the  covariance  does  not  increase  significantly.  However,  Figure  39 
shows  a  graph  of  the  total  position  standard  deviation  as  defined  by 
equation  63.  The  standard  deviation  is  seen  to  rise  during  pass  lb.  Al¬ 
though  this  diverging  tendency  is  not  sufficient  to  cause  filtering  prob¬ 
lems,  the  smoothing  algorithm  as  seen  in  equations  38  through  40  is  ham¬ 
pered  by  this  behavior.  Figure  40  shows  the  velocity  standard  deviation. 
Covariances  generated  by  each  of  the  filters  are  comparable. 

Filter  Effectiveness.  Both  the  WLS  and  the  Kalman  filters  showed  a 
high  degree  of  capability  in  processing  the  selected  data  for  a  polar 
circular  orbit.  Despite  known  limitations  the  models  provided  accurate 
pointing  information.  The  smoothing  algorithm  provided  no  useful  addition 
al  information  during  pass  1  but  functioned  as  expected  during  passes  2 
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Table  XII 

Interpass  F’ointing  Data  Ent  Set  1 


Pass 

Source 

Range 

(meters) 

Range  Rate 
(meters/secotid) 

Elevation 

(deg) 

Azimuth 

(deg) 

la-lb 

WLS 

8S612S 

4784 

33.59 

239.59 

la- lb 

Kalman 

856136 

4790 

33.55 

239.59 

la-lb 

Measure 

859086 

4785 

32.50 

239.52 

lb-2 

WLS 

1726472 

-3848 

10.05 

348.06 

lb-2 

Kalman 

1829049 

-4353 

8.55 

353.  03 

lb- 2 

Measure 

1822414 

-4333 

8.70 

352.48 

2-3 

Kalman 

2338937 

-5107 

2.58 

348.35 

2-3 

Measure 

2334897 

-5099 

2.72 

348.12 
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Chapter  VII 

Actual  Data:  Ent  Set  2 

Data.  Data  for  the  orbit  of  satellite  6100  was  supplied  by  the  1st 
Aerospace  CON  Squadron,  Ent  AFB,  Colorado  for  use  in  this  study.  This 
data  was  processed  at  Ent  using  a  general  perturbations  least  square  error 
processor  (SQEPRI3C).  The  Ent  data  consisted  of  74  observations  of  range, 
range  rate,  azimuth,  and  elevation.  The  results  were  presented  at  an 
epoch  time  corresponding  to  the  time  at  which  the  vehicle  crossed  the 
equitorial  plane  ascending  in  a  northerly  direction.  The  details  of  the 
data  are  presented  in  Table  XIII. 

Radar  Characteristics.  Sensor  344  was  used  exclusively  for  this 
study.  The  parameters  of  this  radar  are  given  by  Table  XIV  as  they  were 
supplied  by  Ent.  The  sigmas  given  are  actually  extremely  high  for  the 
data  of  this  study  and  more  attention  wi-1  be  given  this  problem  later  in 
this  chapter. 

Comparison  of  Orbital  Parameters.  The  orbital  parameters  supplied 
by  Ent  were  of  the  mean  or  general  perturbations  type  (Ref  21:231).  To 
compare  these  parameters  with  those  obtained  by  the  filters,  a  graphical 
analysis  of  the  orbital  parameters  is  presented  in  figures  39  through  48. 
The  results  of  this  analysis  are  presented  in  Table  XV.  It  should  be 
noted  at  this  point  that  the  processor  used  by  Ent  found  the  data  diver¬ 
gent  although  it  obtained  an  orbit  where  the  residuals  were  extremely  small 
and  therefore  the  orbital  parameters  approximately  truly  determined.  The 
parameters  found  by  the  WLS  filter  are  comparable  to  those  obtained  b>  the 
Ent  processor  when  the  divergence  of  the  Ent  processor  is  considered. 
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Table  XIII 

Orbit  Data  Ent  Set  2 
Satellite  #6100  Day  215  Year  72 
=  20  hr  42  min  57.365  sec  (Ref  Appendix  F) 
Start  GMT  =  9  hr  25  min  39.134  sec 


No.  of 

Pass  Obs  Int  Observations 


Duration 

Duration 

Pass 

Gap 

Min  Sec 

Station 

Min  Sec 

i 

j 


1  6  sec  30  3:00  344 

104:34 

2  6  sec  44  4:24  344 
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Table  XIV 

Sensor  #344  Ent  Set  2  Parameters 

Position 

Longitude:  185.89769°W 

Latitude:  52.73267°N 

Height:  86.  meters 

Measurement  Biases 
Range :  +160  m 

Range  Rate:  1  m/sec 
Elevation:  .063° 

Azimuth:  -.019° 

Measurement  Sigmas 
Range :  852  m 

Range  Rate:  12.8  m/sec 
Azimuth:  .426  deg 

Elevation:  .304  deg 
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Table  XV 

Ent  Data  Set  2  General  Perturbations 
Epoch  Time:  Day  215  Hour  10  Min  55  Sec  19.920 


Elements 

Ent-** 

Kalman 

WLS*** 

Semi-Major 

Axis  (E.R.) 

1.178309 

1.1793 

1.1783 

Eccentricity 

.122512 

.1236 

.1245 

Angle  of  Inclination 
(Deg) 

48.4313° 

48.453° 

48.422° 

Line  of  Nodes 
(Deg) 

255.6498° 

255.3691° 

255.3C460 

Arg.  of  Perigee 
(Deg) 

172.2760“ 

171.4976° 

171.9439° 

Period  (Min) 

108.0 

108.20 

108.05 

Height  at  Perigee 
(Km) 

216  Km 

216 

205 

**Failed  Convergence 

***Converged  in  5  iterations 


KALMAN 


KflLMflN 


KflLMRN 


KflLMRN 


Period,  Kalman,  Ent  Set 


WLS  ENT  SET 


ELAPSED  TIIIEtniNUTES) 


WLS  ENT  SET 


WLS  ENT 


ELflPSEO  T I  ME  C MINUTES  ) 

•e  49.  Line  of  Nodes,  WLS,  Ent  Set 
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Table  XVI 


Ent 

Data 

Set  2 

Special  Perturbations 

Epoch  Time: 

Day  215 

Hour 

10  Min  55 

Sec  19.920 

Elements 

Kalman 

Reversed 

Kalman 

Revised 

WLS  Predicted* 

Semi -Major 

Axis  (E.R.) 

1.179962 

1.179216 

1.178903 

Eccentricity- 

.123101 

.>•22853 

.123895 

Angle  of 

Inclination  (Deg) 

48.4691° 

48.4544° 

48.4364° 

Line  of  Nodes  (Deg) 

255.3457° 

255.3036° 

255.2716° 

Arg.  of 

Perigee  (Deg) 

171.5305° 

171.8982° 

171.9950° 

Period  (Min) 

108.294 

108.191 

108.148 

Height  at 

Perigee  (Km) 

221  km 

219  km 

209  km 

Line  of  Nodes 

Rate  (Deg/Day) 

-3.8363 

Arg.  of  Perigee 
Rate  (Deg/Oay) 

3.4751 

Radius  (m) 

8448541 

8432774 

8441018 

Velocity  (m/sec) 

6436 

6444 

6438 

^♦Converged  in  5  iterations 


KflLMRN-  SET 


KflLMRN 


Elapsed 
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However,  the  Kalman  filter  produced  results  which  were  less  accepta¬ 
ble  with  respect  to  the  orbit's  period,  and  angle  of  inclination.  A  major 
cause  of  this  problem  was  found  to  be  attributable  to  the  statistics  supplied 
with  the  radar  data.  While  the  azimuth  sigma  was  given  as  .426°  actual 
analysis  of  the  residuals  produced  by  the  convergent  WLS  filter  showed 
that  a  value  of  about  .025  would  be  much  more  acceptable.  A  similar  dis¬ 
crepancy  of  approximately  a  factor  of  10  was  found  in  the  elevation  res¬ 
iduals.  Repeated  runs  showed  the  WLS  filter  to  be  rather  insensitive  to 
the  values  chosen  for  the  sigmas  with  a  period  of  108.06  minutes  occurring 
repeatably  at  tQ.  Having  revised  the  values  of  sigma  to  correspond  to 
those  found  with  the  residuals,  the  Kalman  filter  program  was  rerun  with 
the  results  presented  in  Table  XVI  and  figures  49  through  53.  The  new 
closer  comparison  with  the  values  obtained  by  the  WLS  filter  are  more 
representative  of  the  true  values. 

Kalman  Smoothing.  The  results  of  the  Kalman  smoothing  of  the  first 
observations  of  each  pass  are  given  in  Table  XVII.  As  was  the  case  with 
real  data,  set  1,  no  significant  results  are  obtained  in  pass  1  since  in¬ 
sufficient  measurements  are  present  to  allow  a  steady  state  value  to  be 
reached.  However,  pass  2  does  achieve  a  steady  state  with  the  values  be¬ 
ing  approximately  those  found  by  reverse  integrating  from  the  final  state 
produced  by  the  Kalman  filter.  Table  XVII  shows  the  smoothing  results. 

Covariance  of  the  Filters.  The  final  values  of  the  covariance  P  are 
shown  in  Table  XVIII  for  the  orginal  Kalman,  the  revised  Kalman,  and  the 
WLS.  No  comparable  values  were  supplied  with  the  Ent  data.  The  covari¬ 
ances  of  the  revised  Kalman  filter  are  plotted  in  figures  54  through  57. 

The  results  are  typical  of  those  achieved  with  both  the  theoretical  and 
Ent  data  set  1. 
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Table  XVII 

Smoothing  Data  at  First  Observations,  Ent  Set  2 
Special  Perturbations 


Kalman 

Reversed 

WLS 

Smoother 

Pass  1  t=0.,  30  observations  used 

for  smoothing. 

Semi-Major 

Axis  (E.R.) 

1.17904 

1.17826 

1.18091 

Eccentricity 

.12385 

.12462 

.12162 

Angle  of  Inclination 
(Deg) 

48.4443 

48.4179 

48.4743 

Line  of  Nodes 
(Deg) 

255.3169 

255.5435 

256.3354 

Arg.  of  Perigee 
(Deg) 

171.6493 

171.7375 

171.152 

Period  (Min) 

108.167 

108.060 

108.425 

Pass  2  t=6454.,  44  observations  used  for  smoothing 

Semi-Major 

Axis  (E.R.) 

1.17934 

1.17828 

1.17940 

Eccentricity 

.12383 

.12458 

.12379 

Angle  of  Inclination 
(Deg) 

48.4513 

48.4184 

18.4505 

Line  of  Nodes 
(Deg) 

255.3304 

255.2564 

256.4090 

Arg.  of  Perigee 
(Deg) 

171.8244 

171.9750 

171.5110 

Period  (Min) 

108.208 

108.063 

108.216 

111 


J. 


GGC/F.E/73-13 


Table  XVI II 

Final  Covariances  Ent  Data  Set  2 


Kalman 

Original 

Kalman 

Revised 

Sigmas 

WLS 

WLS 

Revised 

Sigmas 

p 

XX 

6.23  x  105 

1.05  x  105 

2.88  x  106 

1.00  x  101 

pyy 

2.17  x  106 

5.61  x  104 

1.20  x  106 

9.90  x  10' 

pzz 

7.94  x  106 

1.81  x  105 

6.72  x  106 

1.43  x  10 

p« 

32.45 

1,13 

37.20 

1.15 

p.. 

yy 

47.00 

0.89 

15.83 

0.82 

P. . 
zz 

456.42 

7.01 

386.90 

6.74 

il2 
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Figure  56.  Covariance  of  P^,  Kalman  Revised,  Ent  Set  2 
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Interpass  Pointing  Data.  Due  to  the  small  number  of  measurements 
during  pass  1,  less  accurate  pointing  data  would  be  expected  than  was 
achieved  in  either  the  theoretical  case  or  in  data  set  1.  The  pointing 
data  is  shown  in  Table  XIX  with  an  error  of  about  3°  in  azimuth  and  1° 
in  elevation.  Although  this  would  definitely  limit  the  required  field  of 
scan,  downrange  detection  with  this  pointing  data  will  be  more  difficult. 
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Table  XIX 

Interpass  Pointing  Data  Ent  Set  2 


Pass 

Source 

Range 

(meters) 

Range  Rate 
(meters/secorfd) 

Azimuth 

(Deg) 

Elevation 

(Deg) 

1-2 

WLS 

2789916 

-3650 

183.257° 

24.56° 

1-2 

Kalman 

2792497 

-3655 

183.341° 

24.53° 

1-2 

Measurement 

2712827 

-3492 

180.568° 

25.64° 
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Chapter  VIII 

Conclusions  and  Recommendations 

Conclusions.  Based  on  the  material  presented  in  this  thesis,  as 
well  as  knowledge  gained  throughout  the  study  of  the  near-earth  orbital 
determination  problem  the  following  conclusions  are  drawn: 

I.  The  simulated  orbital  analysis  indicated  that: 

1.  The  Weighted  Least  Square  error  filter  provides  extremely 
accurate  orbital  parameters  in  the  multi-pass  problem  with  long  data  gap 
periods . 

2.  The  Kalman  Filter  also  provides  equally  accurate  orbital 
parameters  in  the  multi-pass  problem  with  minimal  data  and  long  time  gaps 
between  observation  passes. 

3.  The  Kalman  Fixed  Point  Smoother  shows  outstanding  potential 
for  the  determination  of  orbital  parameters  if  adequate  measurements  are 
available  to  cause  the  smoother  to  reach  a  steady  state. 

4.  The  existence  of  both  short  term  and  secular  variations  in 
the  classical  orbital  parameters  were  demonstrated  from  the  application 
of  oblate  earth  modeling. 

5.  No  conditions  were  encountered  which  caused  either  the  WLS 
or  Kalman  filter  to  diverge.  Interpass  covariance  remained  well  within 
the  computational  limits  of  the  CDC  6600.  The  effects  of  the  level  I 
modification  were  not  great  due  to  the  extreme  computational  capabilities 
of  the  CDC  6600. 

6.  The  WLS  filter  took  much  longer  to  achieve  the  same  results 
as  the  Kalman  filter  due  to  the  necessity  to  batch  compile  the  data  and 
repeatedly  integrate  through  that  portion  of  the  orbit  in  which  there  were 
no  measurements.  The  approximate  time  ratio  was  6  to  1. 
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7.  Final  covariances  determined  by  the  two  methods  were  com¬ 
parable. 

8.  Interpass  pointing  data  was  more  accurate  from  the  Kalman 
filter  although  no  great  difference  was  observed. 

9.  Modifications  of  filter  parameters  did  not  improve  the  speed 
or  accuracy  of  the  solution. 

II.  The  actual  data  analysis  indicated  that: 

1.  Both  the  Kalman  filter  and  the  WLS  filter  showed  the  capa¬ 
bility  to  accurately  determine  orbital  parameters  in  a  multi-pass  multi¬ 
radar  environment. 

2.  Kalman  Smoothing  failed  to  demonstrate  the  same  degree  of 
capability  as  in  the  theoretical  case  although  when  steady  states  were 
achieved  the  orbital  parameters  determined  were  accurate. 

3.  The  same  secular  and  short  period  perturbations  seen  in  the 
theoretical  data  appeared  in  the  orbital  parameters  for  the  actual  data. 

4.  Both  polar  circular  and  low  perigee  elliptical  orbits  were 
able  to  be  accurately  determined. 

5.  The  Kalman  filter  was  found  to  be  far  more  sensitive  to 
errors  or  poor  statistics  on  the  noise  covariance  matrix,  R. 

6.  Final  covariances  of  both  filters  were  comparable. 

7.  The  WLS  filter  was  seen  to  diverge  on  both  sets  of  actual 
data  until  a  modification  to  the  process  was  made  whereby  the  first  pass 
data  was  converged  and  then  used  as  the  starting  condition  for  the  multi¬ 
pass  data.  Following  the  modification  neither  data  set  diverged. 

8.  The  Kalman  filter  reached  high  interpass  covariances  on  the 
order  of  10**  during  actual  data  processing.  However,  the  filter  always 
recovered  on  update  and  never  diverged. 
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Recommendations :  The  £ollov;ing  recommendations  are  made  for  con¬ 
tinued  study  in  this  subject  area  and  for  improvement  of  the  existing 
computer  programs: 

1.  Modifications  should  be  made -to  both  programs  to  automate 
timing  decisions  during  integration  routines  in  order  to  maximize  the  use 
of  information  available  from  an  observation. 

2.  Determination  of  0.  should  be  written  as  a  subroutine  and 

So 

added  to  both  programs  to  increase  program  efficiency  and  minimize  re¬ 
quired  input  data.  A  constant  update  of  0  should  be  made  to  achieve 

8o 

maximum  accuracy. 

3.  Both  programs  should  be  modified  to  take  the  universal  data 
format  now  employed  at  ENT. 

4.  The  Kalman  filter  program  should  be  expanded  to  include 
estimation  of  the  ballistic  coefficient,  noise  sigmas,  noise  biases,  and 
net  drag  on  the  vehicle  throughout  the  orbit.  The  estimation  and  applica¬ 
tion  of  the  estimate  of  the  noise  covariance  should  improve  Kalman  filter 
operation. 

5.  A  new  orbital  determination  routine  should  be  written  em¬ 
ploying  the  new  Institute  of  Astronautics  and  Aeronautics  orbital  para¬ 


meters  . 


6.  A  sensitivity  analysis  of  both  the  WLS  and  Kalman  filtering 


methods  of  orbital  determination  should  be  made. 

7.  Additional  actual  data  should  be  applied  to  both  programs 
to  determine  if  the  success  found  in  the  two  cases  tried  in  this  thesis 
can  be  extended. 

8.  Additional  applications  of  the  smoothing  algorithm  should 


be  made  to  determine  if  the  smoother,  when  given  adequate  data,  can 
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perform  to  the  capability  demonstrated  in  the  simulated  orbit  case. 

9.  Both  programs  should  be  analyzed  for  possible  reprogramming 
ici  order  to  reduce  core  requirements  and  decrease  time  requirements. 


( 
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Appendix  A 

Derivation  of  the  Gravitational  Field 
of  the  Earth 

The  equations  which  describe  the  gravitational  field  of  the  earth 
are  both  long  and  complex.  For  this  reason  each  of  the  three  components 
of  the  gravitational  field  expressed  in  the  geocentric  rotating  coordinate 
system  will  be  determined  term  by  term  and  combined  finally  as  a  summation 
of  these  individual  terms. 


Generally,  the  gravitational  potential  is  given  by 


U 


.  2 
ke  m 


1  + 


;(S  ( 

k=2'm=o  '  r^  '  ' 


Ck,mcos(mXE)  + 
Sk,msi"'mXE>)  ) 


(3) 


However  for  this  study  adequate  accuracy  was  considered  achieved  for 

(m), 


U  ■ 


ke2m 


1  + 


til  (!s_ 

k=2'm=o  \  r 


(sin4»)  \  / 

~ J  ^Ck,mcos^  + 

Sk>msin(mXE)  )) 


(65) 


In  order  to  proceed  further  it  is  necessary  to  evaluate  the  associated 

* 

Legendre  Polynomials  from  P2  to  p|.  These  functions  are: 

0  2 
P2  «  ,5(l-3x  j 

P2l  =  3x(1-x2)1/2 

P22  =  3(l-x2) 

P3°  ■  .5x(3-5x2) 

P3J  =  3/2(5x2-1)(1-x2)1/2 

P32  =  15x(l-x2) 


125 


■  '!;■ 


GGC/EE/73-13 


P33  -  lS(l-x2)3/2 

P4°  =  «125x(3-30x2  +  35x4) 

P4J  =  2.5x(1-x2)1/2(7x2  -  1) 

P4"  o  7.5(l-x2) (7x2  -  1) 

P43  =  10Sx(l-x2)3/2 

P44  =  105(l-x2) 2 

P5°  «=  .12S(15-70x2  +  63x4) x 

P5J  =  .125(1-x2)1/2(15-210x2  +  31Sx4) 

2 

P5  *  .12S(l-x2) (-420x  +  1260x3) 
p53  -  .125(1-x2)3/2(-420  +  3780X2) 

P$4  =  . 125(l-x2) 2(7560x) 

p5  ■=  945(l-x2)3^2 

P6  =  •0625C5-105X2  +  315x4  -  231x6) 

=  .0625(1-x2)1/2(210x  -  1260X3  +  1386x5) 
P62  -  .0625(l-x2) (210  -  3780x2  +  6930X4) 

P6^  a  •°625(1-x2)3/2(-7560x  +  27720X3) 

P64  =  . 0625 (1-x2) 2 (-7560  +  83160x2) 

P65  =  10395 (1-x2) S/2x 
Pfi6  =  10395 (1-x2) 3 


(66) 


We  can  then  write 

U  s  _§  ^  1  _  (l~3sin2*)  -  _£3  sin*(3-5sin2*) 

-  _£i  (3-30sin2*  +  35sin4*)  -  _£§.  sin$  (15-70sin2*  +  63sin4*) 

8r  8r5 

-  (5-105sin2*  +  315sin4*  -  231sin6*) 


16r 
6  6 


(67) 


\l2  'Jl  (^-1,n*)  HClc,icost'n>E>  *  SMSin(,"XE>> 
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where  the  summation  expresses  the  longitude  dependent  potential  terms. 
The  remainder  of  this  equation  can  be  expanded  using  the  Legendre  poly¬ 
nomials  just  given.  From  basic  physics  it  can  be  shown  that 

<68> 

By  recalling  that  sin$  =  L  ,  cos$  = 


and 

and 


r  ~  ^x^+y^+z^ 


sin  Xc  =  ■■==- .  ,  cos  X„ 

E  2 
rx^+y^ 


rx^+ yL 


(69) 

(70) 


the  various  required  derivatives  can  be  obtained. 

The  following  trigonometeris  relations  are  also  required: 
sin  2x  =  2sinxcosx 

•y  n  n 

cos  2x  =  cos^x  -  sin^x  +  1-2  sin  x 

t 

sin  3x  =  3  sinx  -  4  sin  x 
cos  3x  -  4  cos3x  -  3  cosx 
cos  4x  =  8  cosS  -  8  coszx  +  1 

7 

sin  4.x  =  (8  cos  x  -  4  cosx)  sinx 

sin  iix  =  5  sinx  -  20  sin3x  +  16  sin5x 

5  1 

cos  5x  =  16  cos  x  -  20  cos  x  +  5  cosx 

C  Z 

sin  6x  =  (32  cos  x  -  32  cos  x  +  6  cosx) sinx 

6  4  2 

cos  6x  =  32  cos  x  -  48  cos  x  +  18  cos  x  -  1 

After  performing  the  required  substitutions,  the  equations  found  in 
SUBROUTINE  GRAVITY  are  obtained. 

Complete  copies  of  the  programs  used  in  this  study  including 
WLS,  the  Kalman  filter,  and  the  data  generator  are  available  from 
AFIT/ENE,  Wright-Patterson  AFB,  Ohio. 
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Appendix  C 

Derivation  of  the  System  Matrix,  F 

The  system  matrix,  F,  is  used  in  both  the  Weighted  Least  Squares 
Filter  and  the  kalman  filter.  As  defined  by  equation  (22) 

F  =  1M.  I  (75; 

3x  1  x=x(k+l;k) 

Now  the  derivative  of  a  matrix  with  respect  to  a  vector  can  be  shown  to 
be  (Ref  18:553).* 


3*1 

3*! 

3*j 

3Xj 

**2 

**6 

3x2 

2k  - 

3x2 

3Xj 

3x2 

3x6 

3X6  3Xfi  _  3X6 

3Xj  3X2  3x6 

Using  the  set  of  state  equations  given  by  equation  (13)  it  can  be  seen 
that 

F14  °  F25  "  F36  “  1 


F.,  =  l_ii  +  uAfi  - 


ie  "  v  *  x^  * 


B  *3DensitX  /2v 


d?  c  - -  v 

42  3x3y 


*  x  *  B  «  3De-niA-*I  /2. 
4  Sy 


c  9ZU  „  „  „  SDensity  ,, 

43  3x3z  4  3z 

F^^  =  -  v-B-Density/2  -  xjj  -B-Density/2v 

F45  *  2“ie  “  B -Density- x4x5/v 
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F46  =  -B- Density x4-x6/v 


F,,  =  -  V  x,  -B-  3Pensity/?, 

3y3x  3x 

F„  *  ♦  U.  2  -  VI  B  8Density  /2. 

52  3y2  ie  5  ay  ' 

F__  =  ^1.  -  v.x,-B-  3Density  /2. 

**  3y3z  3z 

F54  *  "2u^e  -  B*Densityx4*x5/v 

Fgs  =  -vB*Density/2  -  x2j*B-Density/2v 

^56  “  "B'Density-Xg-Xg/v 

F,,  .  ^  -  v-Xc-B.  3Deilsity  /2. 

61  3Z3X  X6  5x 


-  VX, 


,.r.  SDensity  /2> 


F63  =  3^.  -  vxc.R.  3Density  /2. 


F64  =  -B- Dense -x4-x6/v 
Ffi3  =  -B* Dense *x5’Xg/v 

F^  =  -v*B*Dense/2  -  Xg*B*Dense/2v 
where  v  =  velocity 

B  <=  ballistic  coefficient  =  .015 
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Appendix  D 

Measurement  Sensitivity  Matrix  M 

The  measurement  sensitivity  matrix  is  made  up  of  the  partial  deriv 
atives  of  the  measurements  p,  p,  A,  and  E  with  respect  to  the  states  e- 
valuated  at  the  current  best  estimate  of  the  state  of  the  satellite. 


3p 

—  •  • 

3x 

3p 

•  *  sir 

ai .  . 

3x 

li 

•  •  "  r 
32 

3A  # 

3A 

3x 

3z 

3E  > 

3E 

•  •  --  -i 

3x 

3z 

M  =  3hi00 
3xj 


Where  M(l) ,  M(2),  M(3) ,  and  M(4)  are  given  by 

M(l)  =  ((x+X)/p,  (y+Y) / p ,  (z+Z)/p,  0,  0,  o) 

x  p(x+X)  y  _  p  (y+Y)  z_  _  p(z+Z)  x+X  y+Y  z+Z 

P  P  P  ’  p  P  P  ’  P 

M(3)  =  ( l/(p2-z  2)\  ( -x  sin(0)-y  sin(<|>)cos(0) ,  x  cos(0)  - 
\  T  /  \  T  T  T 

yTsin(4>)  sin(<J>) ,  y^osC*),  0,  0,  0  ) 

M(4)  =  ^1/p  -Zy2)1/2')  ^cos($)cos(0)-zT(X-X)/p2,sin(0)cos(<f>) 

»  sin(4»)-z^(z-Z)/p2,  0,  0,  o) 

P 

and  where  the  state  vector  x  =f  x, 1  f  xl 


and  x,  y,  and  z  are  the  negative  of  the  station  coordinates  in  the  geocen 

trie  rotating  frame  and  Xp,  yT,  Zj  are  the  position  coordinates  of  the 
satellite  in  the  tangent  frame* 
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Appendix  E 

Gaussian  Distribution  Curves 

The  additional  graphs  contained  in  this  appendix  are  included  since 
they  were  the  basis  for  the  decision  to  make  25  runs  of  the  Kalman  and 
WLS  programs  and  take  statistical  averages  in  order  to  satisfy  theoretical 
requirements  for  a  Guassian  distribution. 

The  region  of  quantitization  used  on  the  graphs  is  .05,  A  finer  grid 
would  have  led  to  graphs  of  more  erratic  behavior  and  a  less  fine  grid 
would  have  camouflaged  the  true  intent  of  the  curve.  It  is  seen  that  at 
100  points  the  density  appears  to  be  hardly  more  than  uniform.  At  more 
than  500  points,  the  density  approaches  the  familiar  bell-shaped  Gaussian 
density  curve.  However,  it  should  be  emphasized  that  all  of  these  samples 
passed  a  Kolomogorov  Smirnov  test  with  a  significance  level  of  .20  when 
compared  with  a  true  N(0,1)  distribution. 
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Figure  60.  Gaussian  Density  Distribution  100  Points 
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ire  65.  Gaussian  Density  Distribution  700  Points 
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67.  Gaussian  Density  Distribution  900  Points 
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Appendix  F 


Determination  of  0 


The  determination  of  0  is  necessary  so  that  the  inertial  reference 

go 

system  may  be  properly  located  with  respect  to  the  rotating  geocentric 

coordinate  system.  Specifically,  0  is  the  hour  angle  between  the  Green- 

go 

wich  meridian  and  the  vernal  equinox  (inertial  x-axis)  at  0  hrs  GMT  on  the 


date  under  consideration. 


This  angle  would  normally  be  determined  by  use  of  an  ephmeris.  How¬ 
ever,  if  the  hour  angle  desired  is  for  a  date  as  recent  as  yesterday, 
such  an  ephemeris  is  likely  to  be  unavailable.  Therefore  the  following 
approximation  was  used: 

0  =  1.7294468  +  [Days  Since  1950]  *  .017202915  (84) 

°o 

where 


[Days  Since  1950]  =  [yr-50]  *  365  +  INT  (XL-i)  -  12  +  Day  (85) 

Reduce  to  less  than  2u  and  INT  means  take  the  integer  part  of. 

Example 

Find  0  for  190  day  of  1972 
•  8o 

then 

8  «»  1,7294468  +  [(72-50)*365  +  INT  (Zili)  -  12  +  190]* .017202915 

°o  4 

=  1.7294468  +  [22*365  +  17  -  12  +  190]*. 017202915 
=  1.7294418  +  (8225)*. 017202915 
=  1.7294468  +  141.4939759  =  143.2242268 

=  4.99334607  rad  =  286.09765789°  =  19  hr  4  min  23.437884  sec 
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Appendix  G 

Determination  of  Classical  Orbital  Parameters 
Using  the  State  Vector 

The  determination  of  the  classical  orbital  elements 
a,  the  semi-major  axis 
e,  the  eccentricity 
i,  the  angle  on  inclination 
Q,  the  line  of  nodes 
and  to,  the  argument  of  perigee 
from  the  state  vector  composed  of  the  x,y,z  components  of  position  and 
the  x,y,z  components  of  velocity  expressed  in  the  geocentric  rotating 
coordinate  frame,  e,  is  necessary  to  obtain  many  of  the  figures  in  this 
study.  Additionally,  the  parameters 

P,  the  period 

and  hp,  the  height  at  perigee 

are  of  interest  in  many  orbital  problems.  This  appendix  provides  a  method 
for  obtaining  these  orbital  parameters  using  the  state  vector. 

The  first  step  is  to  convert  the  state  vector  to  the  inertial  comp¬ 
utational  frame,  i,  by  use  of  the  directirn  cosine  matrix  found  in  Appen¬ 
dix  B. 

In  the  inertial  frame,  the  radius,  r,  and  the  velocity,  r,  of  the 
vehicle  can  be  expressed  as: 
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where  the  position  states  x,  y,  and  z  have  been  normalized  using  the 
astronomical  distance  unit  r  =  6,378,165  meters  and  the  velocity  states 
have  been  divided  by  the  astronomical  velocity  unit,  vQ  =  7905.376  meters/ 
second. 

Since  r  x  v  »  h  it  is  seen  that 
hx  =  yz-zy 

V 


a  zx-xz 


(88) 


•  • 


\  B  xy-yx 


where  h  is  the  angular  momentum  and  h  ,h  ,h  are  the  x,y,  and  z  components. 

x  y  z 


Also 


h  a  / 


h2  +  h2  +  h2 
x  y  z 


(89) 


since  a  vector  magnitude  is  the  square  root  of  the  sum  of  the  squares  of 
its  components. 

The  eccentricity,  e,  is  then  found  by 


®  =  V2/h2)  ♦  (h/r  -  l.)2 

The  semi-major  axis,  a,  is  given  by 

a  =  h/(l.-e2) (E.R.) 

The  period,  P,  can  be  found  by 

P  a  2ir((a  re)^/io)  *^/60  min  where  a>  =  3.986012x10^ 
The  height  at  perigee,  h,  can  be  found  from 

hp  °  re(a-a*e-l.)  meters 


(90) 


(91) 


(92) 


(93) 


From  geometric  projections  the  following  angles  can  be  determined: 

i  =  cos”*(hz/h)  deg  (94) 

0  =  tan”*(-hx/hy)  deg  (95) 

One  possible  scheme  for  finding  w  is 

to  =  tan”*(z  cos  i/(x  cos  fi  +  y  sin  0))  -  cos’*((h2-wr)/eh^)  deg 

(96) 
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Care  must  be  exercised  on  assigning  the  proper  algebraic  sign  to  the  values 

of  all  of  the  inverse  trigonometric  functions  with  decisions  based  on  the 

sign  of  h_.  (Ref  21:52-53). 
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